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ABSTRACT 


The  conversion  efficiency  of  parametric  amplification 
in  fluids  is  low  because  of  the  low  dispersivity .  A 
discontinuous  change  in  phase  velocity  at  the  boundary  of  a 
waveguide  introduces  dispersion,  which  in  turn  affects 
conversion  efficiency.  It  is  the  purpose  of  this  thesis  to 
develop  from  first  principles  an  analytical  model  which  may 
be  used  to  numerically  predict  the  conversion  efficiency  of 
a  flat-plate,  acoustic  waveguide  given  the  physical 
parameters  of  the  system. 

To  quantify  weak,  finite-amplitude  interactions  in  the 
guide,  the  linear  behavior  of  the  system  is  analyzed  using 
Green's  functions.  Once  the  linear  characteristics  have 
been  determined,  nonlinear  phenomena  are  investigated,  both 
analytically  and  numerically  via  digital  computer  graphics. 
The  physical  parameters  in  the  numerical  examples  are  chosen 
to  correspond  with  materials  used  in  previously  published 
experimental  work  using  cylinders  rather  than  flat  plates. 
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Chapter  I 


INTRODUCTION 


The 

interaction  of 

finite-amplitude 

acoustic  waves 

in 

layered 

media 

differs 

s ignif ica 

ntly 

from  free-field 

behavior  . 

In 

order  to 

inve  s  t iga  t e 

these 

differences. 

the 

pr opert ies  of 

nonl i near 

acoustic 

wave 

interactions 

in 

unbounded 

media  must  first  be  considered.  In  order 

to 

address 

this 

issue,  therefore. 

the  thesis  begins 

by 

comparing 

the 

r e la  t ive 

advantages 

and 

disadvantages 

of 

parametric  and  conventional  linear  transducers. 


1.1  Parametric  Arrays 


According  to  linear  acoustic  theory,  the  propagation  of 
a  disturbance  in  an  isotropic,  homogeneous  medium  is 
governed  by  the  linear  wave  equation: 


(1-1) 


2 

j)4>  «  0 


2 


2  -2  3^ 

where  the  Laplacian  V  and  the  term  Cq  -^2  operate  linearly 

on  the  the  velocity  potential  <p .  Hence,  the  spectrum  of  any 

wave  must  remain  constant  as  the  disturbance  propagates 

through  the  acoustic  medium.  In  order  to  obtain  Equation 

(1-1),  all  second-  and  higher-order  terms  of  the  velocity 

potential  have  been  neglected.  Whether  this  is  a  "safe" 

assumption  or  not  depends  on  the  value  of  the  acoustic  mach 

number  £  and  the  absorption  characteristics  of  the  medium.* 

In  an  inviscid  fluid,  nonlinearity  plays  a  role  at  any  value 

of  £  .  Moreover,  in  a  viscous  fluid,  e=0.1  is  the  limit 

2 

above  which  second-order  theory  begins  to  fail. 

The  governing  finite-amplitude  wave  equation  with 
second-order  terms  included  is  given  in  Chapter  III.  It  is 
sufficient  to  note  here,  that  the  distortion  of  a 
propagating  waveform  and  a  consequent  change  in  its  spectrum 
as  a  function  of  range  are  described  by  these  terms.  For 
example,  an  initially  monotonic,  finite-amplitude 

disturbance  of  frequency  to  ^  in  a  lossless  medium  will 
eventually  become  a  sawtooth  wave,  the  spectrum  of  which 
contains  all  harmonics  of  the  sinusoid.  The  growth  of  these 


1  u  . 
e=- ,  where 

small  signal  speed  of  sound. 


u  is  the  particle  velocity,  and  c  is  the 


2 

A  more  viable  measure  of  the  manifestation  of 
nonlinear  behaviour  in  a  lossy  medium  is  provided  by  the 
acoustic  Reynolds  number  which,  for  plane  waves,  assumes  the 
form  r=8Ek/a  which  gives  a  ratio  of  nonlinearity  to 
absorption  (viscosity)  loss  per  wavelength.  If  T  >1,  then 
nonlinear  effects  become  dominant.  See  also  Rudenko  and 
So luya n ,  p .  II. 


nonlinearly  generated  components  is  at  the  expense  of  the 
amplitude  of  the  fundamental  (energy  is,  of  course, 
conserved  in  a  lnviscid  fluid)  as  shown  in  Figure  1  for  the 
case  of  an  initially  monotonic  wave  in  an  lnviscid 
dispersionless  fluid. 

For  the  case  of  an  initially  bifrequency  waveform 

(frequencies  w  ^  and  not  only  are  harmonics  of  each 

fundamental  generated,  but  interaction  between  the  two 

fundamentals  produces  in t er mod ul a t ion  components  as  well. 

Figure  2  schematically  depicts  the  growth  and  decay  of 

several  components  with  range.  Increased  absorption  at 

higher  frequencies  is  reflected  in  the  sketch  where  the 

lowest  possible  spectral  component,  the  difference  frequency 

oj_  ,  continues  to  propagate  after  the  higher  frequencies  have 

been  absorbed.  Thus,  the  medium  has  the  effect  of  a  low- 
3 

pass  filter. 

The  term  "parametric  array"  refers  to  the  finite- 
amplitude  generation  of  these  in termodula tion  frequency 
components  by  a  bifrequency  source.  The  "length"  of  the 
array  refers  to  the  region  of  interaction  within  which 
energy  is  transferred  from  the  primary  waves  to  the 
nonlinearly  generated  components. 

Parametric  arrays  are  highly  directional.  For  example, 
the  half-power  beamwidth  of  a  nonlinearly  generated 


^This  is  only  true  of  thermo-viscous  fluids. 
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Fenlon,  private  communication,  August  6,  1976 
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difference  frequency  is  considerably  narrower  than  if 
generated  directly  by  a  conventional  piston  projector 
radiating  at  the  sane  frequency;  i.e.,  in  order  to  produce 
the  same  beamwidth  conventionally,  a  much  larger  transducer 
aperture  would  be  required. 

However,  a  major  factor  inhibiting  more  extensive  use 
of  parametric  arrays  is  their  low  conversion  efficiency. 
"Conversion  efficiency"  refers  to  the  ratio  of  energy 
transferred  to  the  difference-frequency  component  versus  the 
amount  of  energy  initially  entering  the  system  via  the 
primary  waves.  In  a  nondispersive  fluid,  transferred  energy 
is  divided  among  all  nonlinearly  generated  components,  only 
a  small  fraction  being  supplied  to  the  difference  frequency. 
However  in  a  dispersive  medium  where  phase  velocity  is  a 
function  of  frequency,  some  components  will  interact 
resonantly,  thus  increasing  in  amplitude  with  range,  while 
others  will  be  excited  asynchronously  producing  the  spatial 
beating  effect  described  in  Chapter  IV.  Under  such 
circumstances,  a  limited  amount  of  energy  will  be 
transmitted  to  those  components  which  are  nonresonantly 
excited,  since  their  amplitudes  never  exceed  a  maximum  value 
where  synchronous  interactions  are  limited  only  by  the 
initial  strength  of  the  primary  fields. 


1.2  Statement  of  Problem 

As  stated  above,  a  dispersive  medium  is  one  in  which 
the  phase  velocity  of  a  monotonic  disturbance  is  a  function 
of  frequency.  Dispersion  may  be  "medium-induced"  (e.g.,  via 
Internal  relaxation  mechanisms  or  inhomogeneities  such  as 
air  bubbles  in  water  or  embedded  in  rubber)  or  "boundary- 
induced",  the  latter  resulting  from  discontinuities  in 
density  and  bulk  speed  of  sound  at  the  interfaces  between 
two  media. 

Figure  3  depicts  a  three-layered  acoustic  medium  where 
the  discontinuities  at  x“a  and  x*=0  represent  the  boundaries 
of  a  flat-plate  waveguide  (medium  II).  Some  of  the  energy 
entering  medium  II  will  be  internally  reflected  at  the  upper 
and  lower  interfaces  as  it  propagates  in  the  positive  z- 
direction.  It  will  be  shown  that  the  speed  at  which  this 
trapped  energy  propagates  through  the  guide  depends  not  only 
upon  which  mode  it  is  in,  but  also  upon  the  frequency  of  the 
disturbance.  Therefore,  if  parametric  interaction  occurs  in 
a  waveguide,  where  dispersion  is  induced  by  the  boundaries, 
certain  spectral  components  will  interact  resonantly  while 
others  will  be  asynchronously  excited. 

It  is  the  purpose  of  this  investigation  to  develop 
expressions  which  may  be  used  to  analytically  determine  the 
extent  to  which  the  boundary-induced  dispersion  of  a 
waveguide  can  enhance  the  conversion  efficiency  of  certain 


II  and  I,  respectively 
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nonllnearly  generated  frequency  components,  namely,  the  sum 
and  difference  frequencies*  The  Investigation  Is  conducted 
In  two  distinct  parts.  First,  exact  analytic  expressions 
are  derived  to  clearly  define  the  dispersive  relationships 
of  an  acoustic,  slow  waveguide  sandwiched  between  two  semi¬ 
infinite  media  of  arbitrary  characteristic  impedances. 
Second,  once  these  relationships  are  established,  the 
effects  of  dispersion  on  specific  nonllnearly  generated 
intermodula tion  components  can  then  be  evaluated  for  various 
physical  systems.  To  make  this  evaluation,  expressions  for 
the  sum  and  difference  frequency  velocity  potentials  of  an 
initially  bifrequency  wave  are  found  via  the  derivation  of 
Green's  functions. 

It  should  be  noted  that  the  equations  representing  the 
velocity  potential  for  the  sum  and  difference  frequencies 
inside  the  guide  (see  Chapter  III)  contain  terms  which  grow 
linearly  without  bounds  as  the  wave  propagates  through  the 
waveguide.  These  components,  which  are  referred  to  as 
"secular  terms",  occur  for  the  case  of  resonant  interaction 
in  a  lossless  medium.  This  physical  inconsistency  Is 
heuristically  explained  in  Chapter  IV.  However,  to  define 
the  bounds  of  resonantly  excited  components,  constrained 
perturbation  theory  may  be  applied  to  the  problem  (e.g.,  the 
method  of  strained  parameters),^  a  task  that  remains  a 
future  research  objective. 


^A.H.  Nayfeh,  Per turba  t ion  Methods  (New  York:  Wiley- 
Interscience ,  1973),  Section  3.1» 


}■ 


1.3  Theoretical  Framework 


The  sketch  shown  in  Figure  3  represents  a  flat-plate, 
acoustic,  slow  waveguide  (medium  II)  of  finite  thickness  "a" 
sandwiched  between  two  semi-infinite,  homogeneous,  fluid 
half-spaces  of  arbitrary  characteristic  impedances  p^Cj  and 
P3C3,  where  p  denotes  density  and  c  denotes  phase  velocity. 
The  fluid-fluid  interfaces  at  x  =  0  and  x=a  extend  infinitely 
in  the  positive  and  negative  y-directions  and  seni- 
infinitely  in  the  positive  z-direction.  The  term  "fluid"  is 
used  in  this  context  to  describe  a  medium  having  a  shear 
modulus  low  enough  to  neglect  the  effect  of  shear  waves, 
thus  permitting  only  comp  res s i ona 1  wave  propagation  to  be 
considered.  Gases,  most  liquids,  and  some  silastic  rubbers 
exhibit  this  property.  It  is  also  assumed  that  the  free- 
field  phase  velocity  in  medium  II,  i.e.,  the  sandwiched 
layer,  is  less  than  that  in  either  medium  I  or  III.  Without 
loss  of  generality,  it  is  assumed  that 


(1-2)  c^c^^ 


where  c^  (i=l,2,3)  represents  phase  velocity  in  medium  "i". 

As  shown  in  Figure  3,  the  wavenumber  vector  k^  normal 
to  the  plane  wavefront  is  decomposed  into  transverse  and 
axial  components  ( and  y^,  respectively)  such  that 


2  .  y  2 

Ki  +  Xi  ’ 


1 
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where  co  represents  angular  frequency  (co=2nf)  •  The 
horizontal  y-component  is  assumed  to  be  zero;  hence,  the 
investigation  is  a  two-dimensional  problem.  The  governing 
Helmholtz  equation  for  a  harmonic  field  exp(-iwt)  in  medium 
II  is 


(1-4) 


(V2  +  k22)  -  0  . 


where  the  Laplacian  V  in  rectangular  coordinates  is  given 
by : 


(1-5) 


v2  .  (*--  0) 

3x2  3z  ’  3y 


and  the  velocity  potential  in  the  frequency  domain  <J>m  is 


(1-6) 


v  =  V<{>^(x, z)  exp(-iu)t) 


or 


(1-7)  v  =  ViJ)  (x,z)  , 

CO  CO 

wh  ere 


(1-8) 


V 


0) 


■'ma; 
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Notice  that  the  unit  vector  in  the  z-direction  k  in  Equation 
(1-8)  is  not  related  to  the  free-space  wavenumber  vector  k 
in  Equations  (1-3)  and  (1-4). 


1.4  Literature  Review 

Nonlinear  wave  interactions  in  a  dispersive  medium  have 
been  investigated  by  Nayfeh  and  Tsai^  in  order  to  evaluate 
the  nonlinear  effects  of  the  gas  and  lining  material  in  a 
duct.  They  determined  a  third-order  uniform  expansion  using 

g 

the  method  of  multiple  scales  to  analyze  the  nonlinear 

effects  on  the  propagation  and  attenuation  of  all  existing 

modes  in  a  two-dimensional,  hard-walled  duct  lined  with  an 

9 

acoustical  material.  Vaidya  and  Wang  utilized  a  second- 
order  expansion  to  determine  the  spectral  energy  transfer  in 
a  lined  hard-  or  soft-walled  duct.  Both  of  the  above 
investigations  were  primarily  concerned  with  nonlinear 
effects  on  attenuation  of  a  fundamental  wave. 


A.H.  Nayfeh  and  M.  Tsai,  "Nonlinear  Acoustic 
Propagation  in  Two-Dimensional  Ducts,"  Journal  o  f  the 
Acoust leal  Society  o  f  Arne  r ic  a  ,  55  (1  974),  1  166-1  1  72  . 

g 

A.H.  Nayfeh,  Perturbation  Methods  (New  York:  Wiley- 
Interscience ,  1973),  Chapter  6. 

9 

P.G.  Vaidya  and  K.S.  Wang,  "Nonlinear  Propagation  of 
Complex  Sound  Fields  in  Rectangular  Ducts,  Part  I:  The 

Self-Excitation  Phenomenon,"  Journal  of  Sound  and  Vibration, 
5  0  (19  7  7)  ,  2  9-4  2  . 
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Acoustic  parametric  amplification  has  been  suggested  by 
Ostrovskii  and  Papilova,*^  who  investigated  the 
amplification  of  a  fundamental  wave  in  an  acoustic  waveguide 
by  injecting  the  second  harmonic  and  utilizing  the 
dispersivity  of  the  guide  to  prevent  higher-order 
interactions.  This  work  was  again  limited  to  rigid  or  free 
boundaries . 

For  linear  sound  propagation,  directivity  enhancement 
of  a  conventional  circular  aperture  transducer  via  an 
acoustic  slow  waveguide  (e.g.,  silicone  rubber  immersed  in 

water)  has  been  observed  experimentally  by  Rogers  and 

11  12 
Trott  and  investigated  numerically  by  King. 

1  3 

Ryder,  Rogers,  and  Jarzynski  experimentally  and 
numerically  investigated  the  radiation  of  a  nonlinearly 
generated  difference  frequency  in  a  silicone  rubber  cylinder 
of  finite  length.  However,  their  results  are  so  completely 


L.A.  Ostrovskii  and  I. A.  Papilova,  "Nonlinear  Mode 
Interaction  and  Parametric  Amplification  in  Acoustic 
Waveguides,"  Soviet  Physics-Acoustics,  19  (1973),  45-50. 

**P.H.  Rogers  and  W.J.  Trott,  "Acoustic  Slow  Waveguide 
Antenna,"  Journal  o  f  the  Acoustical  Society  o  f  America .  56 
(1974) ,  1111-1117. 

1  2 

B.J.  King,  "Numerical  Investigation  of  an  Acoustic 
Slow  Waveguide,"  Journal  o  f  the  Acoustical  Society  of 
America.  62  (1977),  1389-1396. 

1  3 

JiD.  Ryder,  P.H.  Rogers,  and  J.  Jarzynski,  "Radiation 
of  Difference-Frequency  Sound  Generated  by  Nonlinear 
Interaction  in  a  Silicone  Rubber  Cylinder,"  Journal  of  the 
Acoustical  Society  o  f  Amer lea ,  59  (1976),  1077-1086. 
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dependent  on  numerical  analysis  that  no  phenomenological 
predictions  can  be  based  on  it. 


1.5  Outline  of  the  Analytical  Approach 

To  derive  Green's  functions  for  the  problem  under 
consideration,  the  linear  response  of  a  three-layered  medium 
to  a  sinusoidal  point  source  of  unit  strength  in  medium  II 
is  deduced.  The  solutions  thus  obtained  which  are  expressed 
in  wavenumber  space  are  then  transformed  into  frequency 
space  via  complex  integration  techniques.  The  exact 
integral  includes  contributions  from  complex  as  well  as  real 
poles.  However,  as  discussed  in  Chapter  II,  only  real 
poles,  which  represent  trapped  energy,  are  considered.  It 
is  these  discrete  values  of  the  axial  wavenumber  that  define 
the  dispersion  of  the  waveguide  for  guided  modes. 

Distortion  due  to  nonlinear  wave  interactions  in  the 
medium  is  represented  by  second-order  terms  of  the  acoustic 
wave  equation  given  in  Chapter  III.  For  "weak" 
interactions,  these  second-order  terms  combine  to  form  a 
forcing  function  which  may  be  treated  as  a  source 
distribution.  The  analytical  form  of  the  source 
distribution  is  convolved  with  the  Green's  functions  to 
yield  the  sum-  or  difference-frequency  velocity  potential  in 
medium  II  for  guided  modes. 


Finally,  making  use  of  the  derived  solution,  various 
characteristics  of  the  nonlinearly  generated  sum  and 
difference  frequencies  are  investigated  via  computer 
graphics  techniques  as  outlined  in  Chapter  IV. 
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Chapter  II 

GREEN'S  FUNCTIONS 

2.1  Definition 

In  general,  the  Green's  function  represents  the 
solution  of  a  partial  differential  equation  for  a  harmonic 
source  of  unit  strength  satisfying  specified  boundary 
conditions.  As  stated  in  the  introduction,  the  behavior  of 
free,  dilatational  waves  in  an  acoustic  medium  is  governed 
by  the  Helmholtz  equation 

(2-1)  (V2  +  k2)$  =0  , 

0J 

where  <J>  represents  the  velocity  potential  in  the  frequency 
domain.  Consider  a  point  source  of  unit  strength  at  a  point 
(x',z')  in  medium  II  (i.e.,  0<x'<a,  0<z')  oscillating  with 
the  same  frequency  as  above.  The  governing  equation 
describing  the  acoustic  field  at  (x,z)  in  medium  II  is  given 
by 


(V2  +k2)  GjCx.zJx’  ,z')  =  -6(x-x')5(z-z')  , 


(2-2) 


? .  .. 
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where  Gu(x,z|x' ,z')  represents  the  Green's  function  and  the 
Dirac  delta  functions  on  the  right-hand  side  are  defined  to 
be  zero  everywhere  except  at  the  point  x  =  x',  z  =  z'  where  the 
amplitude  becomes  infinite  but  the  magnitude  (or  the  area 
under  the  curve)  is  unity.  Mathematically,  this  can  be 
expressed  by  the  Dirac  delta  function: 


(2-3) 


a)  6  (x)  =0  x  ^  0 


b) 


6  (x)dx  =»  1 


—  CO 


Likewise,  it  can  be  shown  that; 


(2-4) 


CO 

f (x) 5 (x)dx  =  f (0) 


—  CO 


The  Green's  function  G^ ( x  ,  z | x  '  ,  z ' )  can  thus  be  thought  of 
as  the  normalized  response  at  an  observation  point  (x,z)  to 
a  point  source  of  unit  strength  at  (x',z').  It  is  then 
obvious  from  Equation  (2-2)  that  the  Green's  function  will 
be  independent  of  a  particular  source  distribution  since  the 
right-hand  side  of  Equation  (2-2)  represents  a  point  source. 


2.2  Velocity  Potential  in  Terms  of  Green's  Functions 

If  a  volume  source  distribution  Sw(x,z)  at  frequency  ui 
exists  in  medium  II,  the  wave  equation  assumes  the 


inhomogeneous  form 
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(2-5) 


(V2 +k2)f  (x,z)  = -S  (x, z) 

u  u 


From  the  theory  of  linear  differential  operators,  the 

solution  of  Equation  (2-5)  in  the  guide  is  given  by  the 

14 

convolution  of  and  : 


(2-6) 


<t>  (x.z)  = 
0) 


z  a 

'  * 

S  (x’.z')  G  (x,2|x* ,z')dx'dz' 

LJ  U 

0  0 


The  above  solution  is  valid  for  guided  modes  only,  which  are 
discussed  later  in  the  chapter. 

The  power  and  versatility  of  utilizing  Green's 
functions  may  be  seen  in  the  above  equation.  Once 
Gw(x ,z | x ' , z ' )  has  been  determined  for  a  particular  medium 
with  specified  boundary  conditions;  then,  the  response  to 
any  source  distribution  may  be  obtained  via  convolution. 
During  the  remainder  of  this  chapter,  expressions  for 
G w(x , z | x  '  ,  z  '  )  will  be  derived. 

The  first  step  in  the  analysis  is  to  transform 
G  (x , z | x'  , z '  )  into  wavenumber  space.  In  this  manner, 
Equation  (2-2)  becomes  an  ordinary  differential  equation 
rather  than  a  partial  differential  equation.  The  method  of 
sepa r at ion-o f -the-va r iab le s  then  yields  the  assumed  form  for 
the  transformed  Green's  function,  the  coefficients  of  which 
are  obtained  by  matrix  inversion  (or,  in  this  instance,  via 


.  UC.  Lanczos,  Linear  Differential  Operators  (London: 
D.  Van  Nostrand  Co.  Ltd.,  1961),  pp.  206-314. 
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Cramer's  rule).  An  Inverse  transformation  from  wavenumber 
space  to  real  space,  achieved  via  residue  theory,  then 
gives  the  Green's  function  for  medium  II. 

2.3  Derivation  of  G^, 

As  stated  above,  the  Green's  function  satisfies  the 
following  equation: 

H2  2 

(2-7)  ( — «•  +  <  )  G  =  -6(x-x’ )exp(-iXz’ )  , 

dx  X 

where 

00 

’ 

(2-8)  Cx  =  Gx(x|x,,z’)  =  Gaj(x,z|x,,z,)exp(-iXz)  dz  . 

—  00 

Equation  (2-7)  is  a  one-d in ens i ona 1  wave  equation  to 
which  the  following  boundary  conditions  characteristic  of  a 
fluid-fluid  interface  are  applied: 

a)  continuity  of  pressure  at  x=0,a,x' 

b)  continuity  of  normal  velocity  at  x=0,a 

c)  discontinuity  of  normal  velocity  at  x*x'. 

Since  Equation  (2-7)  is  homogeneous  everywhere  except  at 

x=x',  then  G  may  be  treated  as  a  velocity  potential. 

A 

Hence,  the  boundary  conditions  become: 


! 
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3G\  3Gy . 

a)  g-jj-  -  Pj  -g-^— 1  ;  at  x=0,a,x*  for  i,j-l,2,3 

9G  9G 

(2-9)  b)  g^-1  =  ;  at  x*0,a  for  i,j  =  l,2,3 

3G  x'+^ 

c)  limit  _ X  =  -exp(-iXz')  f 

a*  x'-; 

where  p  is  density  in  medium  "i".  The  discontinuity  at  x' 
i 

implied  in  condition  (c)  is  a  consequence  of  the  nature  of  a 
point  source.  It  is  assumed  that  the  point  source  lies  in 
medium  II  (i.e.,  0<x'<a  and  0<z').  Due  to  the  infinite 
extent  of  the  media  along  the  z-axis,  plane  waves  are 
assumed  to  propagate  in  the  positive  z-direction  since  no 
reflected  (i.e.,  backward-going)  waves  are  admitted.  In  the 
transverse  direction,  energy  in  medium  II  is  reflected  at 

the  interfaces  producing  standing  waves  which  are 
represented  as  interfering  plane  waves  propagating  in  the 
positive  and  negative  x-direction.  Energy  escaping  into  the 
semi-infinite  surrounding  media  propagates  away  from  the 
z-axis.  The  appropriate  algebraic  sign  is  chosen  to 
describe  waves  which  decay  to  zero  at  x*»+  00  .  The  free-space 
wavenumber  may  be  decomposed  into  its  horizontal  and 
vertical  components  using  the  following  relationship: 


(2-10)  k± 

where 

K^=transverse  wavenumber  in  medium  "i", 
and  X^=ax^al  wavenumber  in  medium  "i"  . 

Applying  the  method  of  s e para t ion-of -the-var iab les  to 
evaluate  Equation  (2-7)  in  the  various  regions  defined  in 
Figure  3  produces  the  following  solutions: 

(2-11)  a)  x^a  :  G  =A  exp(ie  x) 

A  A 

b)  a2Lxix'  :  G^ = B^exp(i<2x)  + B2exp(-i<2x) 

c)  x'^x^O  :  =  C1exp(ix2x)  +  C2exp(-i<2x) 

d)  O^x  :  =  D2exp  (-iic^x)  , 

where  the  coefficients  A^,  B^,  B2»  C^,  C2>  and  D2  depend  on 

physical  properties  of  the  media.  As  stated  above,  the  sign 

convention  for  the  exponents  in  G  is  chosen  to  avoid 

X 

sources  at  x  =  +  °°  ,  assuming  harmonic  time  dependence 
exp(-iu)t).  It  should  also  be  noted,  that  in  order  to 
preserve  continuity  of  the  wavefront  across  the  boundaries, 
the  axial  wavenumbers,  ,  must  be  the  same  for  all  three 


media: 


(2-12) 


22 


Xl=’  X2  “  X3  *  X 


Applying  the  boundary  conditions  given  in  Equation 
(2-9)  to  the  above  system  yields  a  set  of  six  linear 
equations: 

a)  P^j^exptiKj^a)  =  p2B1exp(ix2a) +p2B2exp(-iic2a) 

b)  p2B1exp(iic2x') +p2B2exp(-iKr2x')  =  p2C1exp(ix2x')  + 

P2C9exp(-iic2x') 


c)  P2C1 + P2^2 = P3°2  * 


(2-13)  d)  K1A1exp(iK1a)  =  <2B1exp(iK2a)  -  <2B2exp(-iic2 


a) 


e)  1*2^  -  iK2C2  *  -iK3D2  »  and 


f)  -iK^expdt^x' )  +  i<2B2exp(-ix2x')  +  iic^exp (ii^x' ) 


i<2C2e3cp  i<2xt  ^  =  exP  (~iXz '  ) 


The  solution  may  be  obtained  via  Cramer's  rule  as  applied  to 
the  matrix  equation  shown  on  the  next  page  [Equation 


(2-14)  ) 
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Expressions  for  the  six  coefficients  Aj,  B^,  C^, 
C^,  and  D2  involve  a  ratio  of  sixth-order  determinants.  A 
computer  program  named  SYMLEQ  (PSU  Computation  Center) 
performs  the  algebraic  manipulations  required  for  the 
solution  of  the  matrix  equation  AX«B,  where  elements  of  the 
matrices  A  and  B  may  be  algebraic  symbols  or  numeric 
constants.  A  discussion  of  the  algorithm  and  general 
guidelines  for  the  use  of  SYMLEQ  are  given  in  a  PSU  program 
guide. ^  Solving  Equation  (2-14)  (via  SYMLEQ)  for  A^,  B^, 
B2’  ^1*  C  2 »  and  °2  an<*  substituting  into  Equation  (2-11) 
leads  to  the  following  expressions  for  Gy  in  each  region  of 
the  transverse  plane: 


Region  I ,  x>a : 


(2-15)  G  =  2exp(-iyz '  )exp  (iic.  (x-a)) 

A  ^ 


[EexpU^x' )  -  Fexp(-ilC2X,  )  ] 


4ic2'i'  Cx) 


H.D.  Knoble,  "Solution  of  Simultaneous  Linear 
Equations  Involving  Matrices  Whose  Elements  are  Symbolic 
Multivariate  (Complex)  Polynomials"  (Program  user's  guide. 
The  Pennsylvania  State  University  Computation  Center,  1971). 
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Region  IIAf  a>x>x': 


(2-16) 


G 

X 


[Aexp(-ix2(x-a))  +Bexp(ilc2(x-a))]exp(ilc2x,)  + 
exp(-ixz'  )  - - 


[ Cexp (-ix2 (x-a) ) + Dexp (i<2 (x-a) ) ] exp(-i<2x ' ) 


Region  IIB,  x'>x>0: 


,,,  {[Aaxp(-ie  (x’-a))  + Bexp (i<_ (x’ -a) ) ] exp(i< _x)  + 

<2-17,  0x  -  e„p (-*,■> - ^ ^ 


[Cexp(-ix2(x'-a)) + Dexp (i<? (x'-a) ) ]exp(-i<2x) y 

4,'nxy 


Region  III,  x<0: 


(2-18) 


Gx  =  2exp(-ixz ' )exp(-i<3x) 


[Gexp (~iic2 (x* -a) )  -  Hexp(i<9 (x'-a) )  ] 
Ak^F(x) 


where  the  following  substitutions  have  been  made  in  order  to 
simplify  the  functional  form  of  the  coefficients: 


,e  l 
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2.4  Inverse  Transformation  of 

The  Green's  functions  in  the  frequency  domain 
G  (x,z|x',z')  are  now  found  by  taking  the  inverse  Fourier 
transform  of  G  : 

A  _____ 

G  (x|x',z')exp(iXz)  dX 

A 

Due  to  the  complex  nature  of  X»  the  integration  extends  over 
the  entire  complex  plane.  The  residue  of  G  at  all  of  its 

A 

poles  must  be  determined  in  order  to  evaluate  Equation 
(2-22)  via  Cauchy's  Residue  Theorem.  However,  as  will  be 
discussed,  the  contributions  from  certain  poles  are 
neglected . 

The  poles  of  G^,  form  a  discrete  spectrum  of  values  of 

the  axial  wavenumber  X»  Real  poles  are  associated  with 

propagating  modes .  It  is  these  modes  which  transport  most 

of  the  energy  through  the  waveguide.  Complex  poles  give 

rise  to  leaky  wave  modes  which  are  more  thoroughly  discussed 

17  18 

by  Kapany  and  Burke,  and  by  Marcuse.  Leaky  modes  are 

t 

t 


C. C.  Ghizoni,  J.M.  Ballantyne  and  C.L.  Tang,  "Theory 
of  Op t ic a 1 -Wav egu id e  Distributed  Feedback  Lasers:  A  Green's 
Functions  Approach,"  IEEE  Jou  rnal  of  Quantum  Electronics,  13 
(1977) ,  843-848. 

^N.S.  Kapany  and  J.J.  Burke,  Optical  Wave  guides .  (New 
York:  Academic  Press,  1972),  pp.  24-34. 

1  8 

D.  Marcuse,  Theory  o f  Dielectric  Optical  Waveguides , 


(2-22) 


(New  York:  Academic  Press,  1974),  pp.  41-46. 
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those  which  continually  radiate  power  from  the  guide  as  they 
propagate  down  the  duct  thus  being  rapidly  attenuated. 
Finally,  in  the  extreme  case,  imaginary  poles  represent 
eva  nes  c  ent  modes  whose  amplitudes  decrease  exponentially 
with  range.  The  effects  of  evanescent  modes  and  leaky  wave 
modes  are  thus  confined  to  ranges  very  near  the  source. 
Hence,  considering  only  propagating  modes  (or  those  for 
which  X  is  real)  and  ignoring  nearfield  effects  permits  the 
integral  (Equation  (2-22)]  to  be  approximated  as  a  sum  of 
residues  of  poles  on  the  real  axes  only. 

The  real,  axial  wavenumber  X  may  range  anywhere  from 
zero  to  infinity  for  each  frequency.  Since  it  is  assumed 
that  C2>Cj>C2»  then  the  free-space  wavenumbers  are  related 
by  k^<kj<k2«  Also,  from  Equation  (2-10),  the  magnitude  of 
the  transverse  wavenumber  in  medium  'i'  is  given  by: 


W^2-*2 


;  i =  1,2,3 


For  X  in  the  range  0<X<k^,  all  transverse  wavenumbers 
are  real  and  it  can  readily  be  seen  from  Equation  (2-11) 
that  this  represents  a  radiation  mode.  For  y  in  the  range 
k^X^k^.  the  transverse  wavenumbers  Kj  and  <2  are  real, 
whereas  is  pure  imaginary.  Hence,  this  mode,  which  is 

confined  to  media  I  and  II  with  energy  decaying 
exponentially  in  medium  III,  is  therefore  called  a  substrate 
mode.  For  k,<X<k2»  as  seen  from  Equation  (2-23),  and  k., 


are  imaginary  and  <2  is  real  so  that  all  energy  propagating 
under  these  conditions  is  restricted  to  medium  II,  hence, 
the  term  guided  mode.  Propagation  of  guided  modes  is 
characterized  by  total  internal  reflection  at  both  upper  and 
lower  interfaces.  Finally,  for  x>^2»  all  transverse 
wavenumbers  are  imaginary  so  that  no  propagating  wave  can 
exist.  The  term  evanescent  mode  is  used  to  describe  this 
condition  under  which  all  energy  decays  exponentially  in  the 
direction  of  wave  propagation.  Regions  for  each  mode  type 
are  depicted  in  Figure  4  which  shows  the  relationship 
between  angular  frequency  0)  and  axial  wavenumber  X* 


2.5  Discussion  of  Poles  of  G ^ 


In  order  to  evaluate  the  integral  defined  by  Equation 
(2-22),  it  is  necessary  to  locate  the  poles  of  the  integrand 
and  evaluate  the  residues.  Since  G y~  F (x)/^(x)  the  roots  of 
the  equation  4'(x)=0  define  the  poles  of  G  .  Setting  the 

A 

function  4'(x)=0  results  in  the  following  transcendental 
equation: 


(2-24) 


tan (k  ^a)  = 


-i(p2P3x1x2  +  P1P2k:2<3^ 


^P1P3‘2 


+  P2  <i<3 


where  the  transverse  wavenumbers  are  given  in  Equation 
(2-23).  Equation  (2-24)  has  two  possible  solutions  for  real 
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X:  either  k  and  are  all  pure  imaginary,  or  <2  is 
real  and  <  and  tc^  are  imaginary.  As  discussed  earlier,  the 
former  case  does  not  permit  a  propagating  wave  to  exist. 
Hence,  for  real  X»  and  <3  are  imaginary  with  <2  real; 
that  is,  X  is  in  the  range  k1<X<k2*  Under  these  conditions, 
only  guided  modes  may  propagate  as  previously  discussed  and 
shown  in  Figure  4.  The  transverse  wavenumbers  may  be 
redefined  to  satisfy  the  above  conditions: 


(2-25) 


K=< 

2  2m 


K-  =  IK 

3  3ra 


where  the  subscript  "m"  denotes  mode  number  as  well  as 
indicating  a  positive,  real  quantity.  Depending  on  the 
magnitude  of  the  argument  of  the  tangent  function 
(i.e.,  <2a ) »  Equation  (2-24)  has  "m"  roots  for  any  given 
frequency,  where  the  mth  pole,  +Xm  ,  satisfies  the  following 
transcendental  equation  obtained  by  substituting  Equation 
(2-25)  into  Equation  (2-24): 


q2P3K1mK2m 


PlP3K2m 


+  PlP2*2mK3m 
'P22"lm<3m 


(2-26) 


tan(K2ma)  = 
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It  is  possible,  when  dealing  with  a  ratio  of  functions,  for 
a  zero  to  cancel  a  pole  and  thus  eliminate  the  apparent 
singularity.  This  is  the  case  for  one  root  of  Equation 
(2-26);  in  particular  Xm“k2*  is  discussed  in  Appendix 
B  (see  Plane-Wave  Mode). 

Now,  the  integrand  of  Equation  (2-22)  is  more  simply 
expressed  in  the  following  form: 


(2-27) 


G  (x,z  x'  ,z')  =r- 
u  m 


G  (x:x' ,z')exp(iXz)  dy 

A 


exp(i(z-z*)x)dx 


where  the  function  F (x )  may  be  deduced  from  expressions  for 
G^  given  in  Equations  (  2  - 1 5  )  -  (  2  -  1  8 )  and  M*(x)  is  defined  in 
Equation  (2-21). 


2.6  Contour  Integration 

As  previously  stated,  the  poles  of  interest  lie  on  the 
real  axis.  By  convention,  the  time  factor  has  been  chosen 
to  be  exp(-iu)t).  Therefore,  real  poles  greater  than  zero 


correspond  to  waves  propagating  in  the  positive  z-direction. 
All  physical  problems  involve  at  least  a  small  amount  of 
damping  which  is  accounted  for  by  including  an  imaginary 
term,  of  appropriate  algebraic  sign,  in  the  axial 
wavenumber.  It  is  apparent  that  for  waves  propagating  in 
the  positive  z-direction  (i.e.,  for  z>z'),  the  imaginary 
part  of  must  be  greater  than  zero  to  satisfy  the 

condition  that,  as  time  increases,  the  amplitude  will 
decrease  to  zero  as  z  -->  <» .  Similarly,  for  z<z',  the 
imaginary  part  must  be  less  than  zero.  The  inclusion  of 
damping,  therefore,  shifts  the  positive  poles  into  the  upper 
half,  and  the  negative  poles  into  the  lower  half  of  the 
complex  X  plane . 

A  section  of  the  contour  of  integration  for  Equation 

(2-27)  is  depicted  in  Figure  5.  Branch  cuts  must  be  made 

/ ~2  2 

due  to  the  multivalued  functions  involved  (i.e.,  X=  *k2  ~K2  )• 
Contributions  to  the  integral  along  these  cuts  have  been 
neglected  for  the  current  investigation.  To  obtain 
exponential  decay  of  the  semicircular  contribution  to 
Equation  (2-27),  the  following  must  be  true: 

(2-28)  limit  |  exp(i(z-z ’ )x)  |  = 0 

!x|  +0 

This  says  that  the  path  of  integration  in  the  upper  half 
plane  applies  to  z  >  z  '  and  includes  the  poles  X  greater  than 


zero.  Then,  by  Cauchy's  Residue  Theorem,  the  Green's 
function  [Equation  (2-27)]  for  forward-guided  modes  is 


given  by  a  sum  of  residues  at  the  positive  real  poles  X  : 

m 

00 


(2-29) 


GJx.zjx'.z*) 


_1_ 

2tt 


Gxexp(ixz)dx 

— oo 


2iri  I  RES  [— ^ 


G  exp(ixz) 


m=l 


2ir 


■;  x  -  x>  o] 

m 


=  I  G 
m=l  “ 


where  G+  represents  the  Green's  function  for  the  mth 
m 

forward-guided  mode. 


2.7  Evaluation  of  Residues 


The  residue  of  a  quotient  of  functions  which  has  a 
first  order  pole  at  X")^  is  given  by: 


(2-30) 


;  X -Xm]- limit  [f^J 


m 


x-*-  X. 


m 


where  V  (X  )  has  been  differentiated  with  respect  to  X  • 
Expressions  for  G+  are  obtained  via  Equations  (2-30)  and 


* 
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(2-29)  together  with  Equations  ( 2- 1 5 ) - ( 2-2 2 )  which 

originally  defined  Gy. 


Region  1,  x>a: 


(2-31)  G  =2iexp(ix  (z-z'))exp(-x  (x-a)) 

m  m  im 


[E  exp(ix  x')  -  F  exp(-ic,  x') ] 
m  Zm  m  2m 


W 


Region  I IA ,  a>x>x 
(2-32) 


{[A  exp(-ix_  (x-a)) 
CT  -  iexp (ix  (z-z'»  — 


w 


B^exp ( in 2m (x-a) ) ] exp ( ix ,bX ' )  +  [ Cmexp (-ix  2m (x-a) ) 


2m  m 

W 


'2m' 


Dmexp(ix2m(x-a)) ]exp(-ix7mx’ ) 


Region  IIB,  x'>x>0 

{ [A  exp(-ix  (x’-a) ) 

(2-33)  Gm  =  iexp(ixm(Z-z’))  - 4<“r(7") - + 

Zm  m 


Bm'^(lr2m(x'‘*))laxP(i,'2m1')+ICmexI,(-i,c2m<K,-*>) 


D[nexp(iic2ra(x,-a))]exp(-iK2nix) } 

V  <*m> 


Region  III,  x<0 


[G  exp(-ix2m(x’-a)) 

G  -2iexp(iX  (z-z’))exp(x  x) - 

tn  *  m  jin 


4x2»'”<Xm> 


W 


i 


(2-34) 
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where  the  following  expressions  are  obtained  by  substituting 


Equation 

(2-25) 

into  Equations 

(2-19)  and  (2-20) 

Le  t 

(2-35) 

a  = 
m 

~°2  KlmK3m 

.  3 

e  =  xp.  K, 

m  2  2m  3m 

b  ■ 
m 

iplP22K2mK3m 

-  2  2 
fm=P2  P3K2n 

c  = 
m 

ip22p3KlmK2m 

gm”ip23<lmK2n. 

d  - 
m 

Plp2P3K2m2 

h  -  p ,pA,  2 

m  1  2  2m 

Then 

) 

A  = 
m 

-a  +b  +c  -d 
m  m  m  m 

and  E  =  e  - f 
id  mm 

(2-36) 

B  = 
m 

a  +b  -c  -d 
m  m  m  m 

F  =  e  +f 
m  mm 

C  = 
m 

a  — b  +c  -d 
m  m  m  m 

G  *  g  -h 
m  m  m 

D  = 
m 

-a  -b  -c  -d 
m  m  m  m 

H  =  g  +h  , 
m  °m  m 

wh  ere 


(2-37) 


'lm 


'2m 


■ \J *22-X» 

■\A 


2  2 

K3m=\/Xm  "k3 


The  expression  for  (X  )  is  derived  in  Appendix  A,  the 
results  of  which  are  restated  here: 


(±xm)[(P2^ 


— — - +  2p  p-p  +p  p  k  a  + 

Klm*3m  123  2  3  1® 


2  2  k2 

plp2  <3ma)sin(K2ma)  +  (plP2P3K2ma  +  p2  P37“ 

lm 


(p22p3Klm  +  P23Klm’C3ma  +  Plp22K3m) 


2  K2mx 


'2m 


+  P1P2  77)cos(K2ma)] 
3m 


By  inspection,  it  can  be  seen  that,  for  x  =  x' ,  in  region 

IIA  is  equal  to  G+  in  region  IIB. 


2.8  Boundary-Induced  Dispersion 

It  is  evident  from  the  form  of  Equation  (2-29)  that  the 
Green's  function  is  expressed  as  an  eigenfunction  expansion 
of  guided  modes  where  the  poles  of  Gy(x|x',z'),  namely  X  , 
are  also  eigenvalues  of  the  system.  Numerical  methods  used 
to  evaluate  these  eigenvalues  will  be  discussed  in  Chapter 
IV. 

Once  the  physical  parameters  of  the  system  are  known 
(i.e.,  density  and  bulk  phase  velocity  of  media  I,  II,  and 
III  and  thickness  of  medium  II),  the  variation  of  axial 
wavenumber  Xffi  with  frequency  for  each  mode  may  be 
determined.  For  example.  Figure  6  shows  frequency  as  a 
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function  of  axial  wavenumber  for  the  first  few  modes  of  an 
infinite  flat-plate  of  silicone  rubber  10  centimeters  thick 
immersed  in  water. 

Several  aspects  of  this  graph  will  be  discussed  later 
in  Chapter  IV.  However,  one  important  point  should  be  made 
at  this  time.  In  a  no nd i spe r s ive  medium,  that  is  one  in 
which  phase  velocity  is  independent  of  frequency,  wavenumber 
and  frequency  are  linearly  related.  Figure  6  illustrates 
the  fact  that  for  any  given  mode,  low-frequency  disturbances 
tend  to  propagate  at  the  higher  phase  velocity  of  the 
external  media  (i.e.,  I  and  III).  Contrariwise,  the  speed 
of  higher  frequency  components  asymptotically  approaches 
that  of  the  waveguide  (medium  II).  This  phenomenon  is 
referred  to  as  "boundary-induced  dispersion"  and  will  be 
discussed  in  greater  detail  in  Chapter  IV. 


FREQ  CKHZ} 

0.00  5.00  10.00  15.00  20.00 


MEDIUM  PHASE  VEL.  DENSITY  THICKNESS 


MODE  CUTOFF (kHz) 


®  1  0.00 

0.10  2  6.70 

00  3  13.41 


*SI  UNITS 


I  1500.0  1000.0 

II  1000.0  1000.0 

III  1500.0  1000.0 


m=3  2  1 


AXIAL  WAVENUMBER  Cl/M} 


Figure  6.  Dispersion  relationship  for  first  few  guided  modes  in  a 
0.1  m. -thick  slab  of  silicone  rubber  immersed  in  water. 
Isospeed  lines  at  X  =  k^,  k^  represent  bounds  for 
trapped  energy.  Dashed  lines  indicate  cutoff  frequencies 
for  modes  m=2,3. 


M  u  I  ■  >  (Iijni 
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Chapter  III 


NONLINEAR  WAVE  INTERACTIONS  IN  LAYERED  MEDIA 


3. 1  Introduction 


The  homogeneous  wave  Equation  (2-1)  is  a  linear, 
partial  differential  equation  obtained  by  neglecting  all 
terms  of  second-  or  higher-order.  If,  however,  convective 
and  medium  nonlinearities  are  taken  into  consideration, 
inclusion  of  second-order  terms  leads  to  the  following 
nonlinear  wave  equation,  where  viscoelastic  and  other  loss 
mechanisms  have  been  neglected: 


(3-1) 


v  2 

3t 


Co2v2)*  =  fi[!v,!2+l^2(^)  ] 

2co 


wi  th 

B  19 

y=ratio  of  specific  heat  in  gases  (1+-  1°  liquids) 

<J>  =the  velocity  potential. 


^R.T.  Beyer,  "Parameter  of  Nonlinearity  in  Fluids," 
Journal  of  the  Acoustical  Society  o  f  America .  9  (  19  60), 
719-721. 
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Equation  (3-1)  is  the  second-order,  nonlinear  wave 
equation.  The  source  term,  or  right-hand  side  of  this 
equation,  contains  the  square  of  first-order  acoustic 
fields.  Hence,  an  originally  monotonic  sound  field  will 
result  in  a  second-harmonic  source  term  over  distances  which 
are  considerably  smaller  than  the  critical  range  at  which 
shock  formation  occurs.  The  new  field,  no  longer  being 
monotonic,  is  again  squared  over  the  next  incremental 
distance  giving  rise  to  interaction  between  the  fundamental 
and  second  harmonic,  which  in  turn  produces  a  third-harmonic 
field  component.  In  an  inviscid  medium,  this  continuous 
process  eventually  leads  to  the  formation  of  a  shock  wave, 
the  spectrum  of  which  contains  all  harmonics  of  the 
fundamental . 

Similarly,  bifrequency  excitation  leads  not  only  to  the 
generation  of  harmonics  of  each  primary  frequency,  but  also 
to  intermodulation  frequency  components  which  result  from 
interaction  between  the  primary  fields.  Hence,  the 
principle  of  linear  superposition  no  longer  holds.  Of 
particular  interest  are  the  sum  and  difference  frequencies 
("denoted  by  >  f°r  a’i  >  ^  )  * 

In  general,  various  loss  mechanisms  inhibit  the 
formation  of  weak  shock  waves  in  viscoelastic  fluids.  Such 
loss  mechanisms  are  included  via  complex  wavenumbers  in  the 
current  investigation.  It  is  assumed  that  the  nonlinear 
interaction  is  weak  enough  to  allow  the  velocity  potential 


•  «-•  ^ 


on  the  right-hand  side  of  Equation  (3-1)  to  be  approximated 

by  a  linear  superposition  of  the  primary  fields.  This 

approximation,  commonly  made  in  weak  finite-amplitude 

20 

acoustics,  has  been  more  thoroughly  discussed  by  Fenlon. 
Using  the  following  Fourier  transform  relationship, 


(3-2) 


where 


(3-3) 


4>(x,z,t)  =  <J>  (x,z)exp(-iiot)dw 

U)  * 


(x,z)  =  -T—  <)>(x,z,  t)exp(iut)dt  f 


the  primary  fields  can  then  be  expressed  via  the  eigenmode 
expansion  at  the  two  frequencies  ^  and 


(3-4) 


$  (x,z)  =  £  exp(i(xUlz  +  ►.-“^-x) ) 


(3-5) 


q  q 


<f>  (x,z)  =  £  Ca’2exp(i(xaJ2z +  icW2x.)) 

“2  s=— 00  s  s  q 


Thus  for  weak,  bifrequency  wave  interactions  in  medium  II, 
if  Equation  (3-2)  is  substituted  into  Equation  (3-1),  and 
only  those  terms  involved  in  sum-  or  difference-frequency 


F.H.  Fenlon,  "On  the  Performance  of  a  Dual  Frequency 
Parametric  Source  Via  Matched  Asymptotic  Solutions  of 
Burg-ers'  Equation,"  Journal  o  f  the  Acoustical  Society  o  f 
America,  55  (1974),  35-46. 


..-'.-VXW:  A.  ./v  -* 


generation  are  retained,  the  resulting  inhomogeneous 
Helmholtz  equation  at  these  frequencies  is  given  by: 


(3-6) 


(V2  +  k+)4> 
Z 


iU+  .  (*>  V_T 

l7*u  •  7*o»  +  klVm.V  ] 

1  2  *  ^ 


where  4>  represents  the  "linearized"  velocity  potential  in 
to . 

medium  II,  and  k.  represents  the  free-space  wavenumbers  in 
medium  II  for  primary  frequency  (ii± ,  W±  being  the  sum  or 
difference  frequency  (i.e.,  w±  =  w2 )  .  2  1  The  superscript 
"*"  in  parentheses  which  appears  in  the  velocity  potential 
in  Equation  (3-6)  implies  that  the  complex  conjugate  only 
applies  to  the  difference  frequency. 

That  is , 


♦<*> 

=  <f> 

u>2 

for 

S  (x ,  z) 
w+ 

s 

* 

=  4> 

W2 

for 

S  (x,z)  . 
co_ 

3.2  Frequency-Domain  Solution 

The  second-order,  nonlinear  wave  equation  in  the 
frequency  domain  is  given  by  Equation  (3-6).  It  may  be 
seen,  by  comparing  this  equation  with  the  general  form  of 


21f.H.  Fenlon,  private  communication,  October  17,  1979 


45 


the  inhomogeneous  Helmholtz  equation  given  in  Equation 
(2-5),  that  the  nonlinear  terms  which  are  grouped  on  the 
right-hand  side  of  the  equation  can  be  thought  of  as  a 
forcing  function,  or  volume  source  distribution  in  the 
guide.  This  being  the  case,  the  velocity  potential  in 
medium  II  for  a  particular  frequency  is  given  by  Equation 
(2-6)  : 


(2-6) 


^Oc.z)  = 


z  a 

r  t 


i  j 
0  0 


s  (x'.z') 

Ui 


G  (x ,  z  j  x 1 ,  z 1 ) dx 1 dz 1 


Hence,  in  order  to  evaluate  Equation  (2-6)  for  the  sum-  or 
difference-frequency  component,  the  source  distribution 
(x,z)  must  be  determined.  By  inspection  of  Equation 
(3-6),  the  latter  can  be  expressed  as 


(3-7) 


S  (x,z) 

Ul 


-  (*)  Y_!  (*) 

7*  +  (V^)  k  k  *  4  ] 

“?  ~  1  2 


[see  Equation  (3-6)  for  definition  of  terms). 

The  first  step  in  determining  <j>  therefore,  is  to 

W± 

evaluate  the  source  distribution  in  Equation  (3-7).  The 
next  is  to  convolve  the  resulting  function  with  the 
appropriate  Green's  function  and  perform  the  integration 
over  x'  and  z'  for  a  typical  mode. 


Td*4*wl*  • 
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3.3  Source  Distribution 


The  sum-  or  difference-frequency  source  distribution 

for  weak  finite-amplitude  interaction  in  medium  II  is  given 

2  2 

in  Equation  (3-7).  The  primary  waves  are  represented  as 
eigenmode  expansions  (see  Equation  (3-6)]  and  hence,  the 
scalar  product  of  the  primary  velocity  vectors  is  given  by 
the  double  summation: 


(3-8) 


(*)  »  (*) 

V<6  •  Vi  =  +  Z  Z  C  1C  1  [c  ^ 

V(D„  q  s  q  s 

1  2  q>  s — " 


ul*-t*>2 


Id.  “9, 

+  \W] 


0),  id-  .  Id.  , , 

exp{i[(x.1  ±  X_2)z  +  1  ±  <s 

q  s  4  3 

Likewise,  the  product  of  the  primary  wave  velocity 
potentials  is  simply 


(*) 


<d.  ld_ 


(3—9) 


®  td1  td^(*)  Id.  Id-  w- 

=  Z  Z  C  lC  exp{i[(x„  +Xc  )*+(«  +<  )*]> 

W1  w2  q,s—  q  S  q  "  q  "  S 


Substituting  Equations  (3-8)  and  (3-9)  into  Equation  (3-7) 
thus  gives  the  following  form  for  the  sum-  and  difference- 
frequency  source  distribution  at  the  source  point  (x',z#): 


2  2 

It  is  assumed  that  all  nonlinear  wave  interactions 
occur  within  the  guide  (i.e»,  0<x<a  and  z>0);  therefore,  the 
velocity  potentials  and  wavenumbers  for  the  primary  fields 
correspond  to  medium  II.  With  this  in  mind,  the  notation 
becomes  slightly  less  laborious. 


(3-10) 
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iu. 


S  (x',z') 


U) 


— ?  E  I  C 
2  q 

Cj  q,s=-“ 


“i  wo(*)  _  wi  “i  “  , 

•  V2  Wx.V!*<; Kj !)  +  (^i: 


■q  s 


q  s 


-)k.k 
2  1 


where 


w,  u-  ll).  ti)0 

exp{i[(x  +xc  )*'  +  (K  ' 

M  S  CJ  S 


m  =  CO  ±  0) 

±  12 


ki  =a)i/c2  ( 1 »  2 ) 


K^'transverse  wavenumber  for  mode  '  q'  , 

frequency  oj  ,  medium  II,  m=l,2. 
m 

va'm  =axial  wavenumber  for  mode  'q', 

Aq 

frequency  (u  ,  medium  II,  m=l,2. 

Cq™ =we ight ing  coefficient  for  mode  ' q' , 

frequency  ;  determined  by  actual  source 

distribution  at  face  of  waveguide  (i.e.,  z=0) 

B  2  3 

y=ratio  of  specific  heat  in  gases  ( 1  +  -  in  liquids). 


Equation  (3-10)  defines  the  distribution  of  source 

points  which  contribute  to  the  growth  of  the  nonlinearly 

generated  sum-  and  d i f f e r enc e- f r eque nc y  fields  inside  the 

waveguide.  These  points  are  referred  to  as  "virtual" 

sources  because  they  act  as  point  sources  scattered 

throughout  the  region  of  interaction  generating  the  harmonic 

24 

and  intermodulation  frequency  components. 

2  3 

R.T.  Beyer,  "Parameter  of  Nonlinearity  in  Fluids," 
Journal  o  f  the  Acoustical  Society  o  f  Amer lea .  9  (1960), 

719-721. 

2  4 

See  e.g«,  H.O.  Berktay  and  C.A.  Al-Temimi,  "Virtual 
Arrays  for  Underwater  Reception,"  Journal  of  Sound  and 
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3.4  Evaluation  of 

Recall  that  the  Green's  functions  for  forward-guided 
modes  in  each  region  are  given  by  a  sum  of  residues  at  real 
poles,  which  is  analogous  to  an  eigenmode  expansion: 


(2-29)  G  (x,z|x',z')  =■  £  G**" 

^  m=  1  m 

Expressions  for  G*  are  given  in  Chapter  II  (Equations 
( 2-3 1 )- (2-34 )] •  Multiplying  the  Green's  functions  for 
region  I IA  [Equation  (2-32)]  by  the  source  distribution 
given  by  Equation  (3-10)  thus  leads  to  the  following 
expression  for  the  velocity  potential  in  medium  II  as  per 
Equation  ( 2-6 ) : 


(3-11) 


4  (x,z) 

u+ 


w,  “  ”  U),  CO  (*) 

-  (_=_)  Z  Z  C  J  R  X  Z 
4c  ^  m=i  q ,  s=-°°  q 


where 


a-1*.  .  -  .  wl  w2  ,  “l  w2. 

R  =  (~)klk2  +  (*q  Xs  K  > 


•fI±T’(xI±) 


V 


<  —  =  transverse  wavenumber  in  medium  II  (see  text) 
at  frequency  u.  for  mode  ra 


and 


Vibration,  9  (1969),  295-307. 
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(3-12) 
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X  ■  transverse  component 

U4.  "+ 

“  IAmexP(_ilcm~(x_a))  +  3_exp(iic_—  (x-a))] 


m 


“4 

m 


exp(iA+x')dx' 


■  w+  (o 

+  [Cmexp(-i»cm—  (x-a))  +  D^expd^—  (x-a))  ] 


m  m 


exp(iAxx')dx'  , 


(3-13) 


where 


Z  =>  axial  component 
z 


“+ 

=  exp(rx  -  ) 

ID 


exp(iAzz' )dz' 


.+  W1 
A  =  K 

+ 

“2 

K 

+ 

<  — 

X 

q 

S 

m 

-  “l 

A  =  K 

+ 

K 

<  — 

X 

q 

s 

m 

A 

A  =  v 

+ 

“2 

X 

X  “ 

z 

q 

s 

Am 

and  A  ,  B  ,  C  ,  D  are  all  defined  in  Equation  (2-36), 
m  in  in  id 

T' (X^1)  being  given  by  Equation  (2-38). 

The  velocity  potential  is  thus  represented  as  a  triple 
summation  over  ra,  q,  and  s,  where  m  represents  the  mode  in 
which  the  sum  or  difference  frequency  propagates,  and  q  and 
s  are  the  modes  excited  by  the  primary  field  sources  tOj  and 
oj2  at  z  =0 . 

In  order  to  examine  the  contribution  of  a  particular 
set  of  modes,  terms  for  positive  and  negative  q  and  s  should 
be  included  for  each  integer  value  of  m,  since  waves  travel 
in  both  the  positive  and  negative  x-directions  Inside  the 


guide.  Therefore,  for  fixed  values  of  m,  q,  and  s,  the  x- 
component  of  Equation  (3-11)  involves  two  terms  designated 
by  a  subscript  (i.e.,  X_^_  and  X _ )  to  distinguish  between  q 
and  s  being  greater  or  less  than  zero.  Then  since 


(3-14) 


b 

exp(iar)dr  =>b  exp(iab/2)sinc(ab/2) 

0 


where  sine (x) =sin (x) /x ,  the  transverse  and  axial  components 
of  Equation  (3-11),  i.e..  Equations  (3-12)  and  (3-13), 
become : 


(3-15)  X+  =  ( J+K)  exp  ( iKm~x)  +  (L+M)  exp(-ix“ 

where 

3=  B  a  sine  (A  a/  2)  exp  (— i  Ctc  ^a-A+a/2)) 
m  x  m  x 

w+ 

K=  D  a  sine  (A  a/2)exp(-i(ic  -a  -  A  a/2)) 
m  x  m  x 

L=  A  a  s inc ( A  a/ 2 ) exp ( i (<  :!:a  +  A+a/2)) 
ra  x  m  x 

M=  C  a  sinc(A  a/2)exp(i(<  — a + A  a/2)) 

m  x  m  x 

and 


...  •  .'v-:  ‘•♦-vtr  _  .  -.-tt.- 
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“+  “+ 

(3-16)  X_=  (U+V)exp(i*  — x)  +  (W+Y)exp(-iic  — x)  , 

in  m 

where 

“+ 

U=  B  a  sinc(A  a/2)exp(-i(*  — a + A  a/2)) 
m  x  m  x 

+  . 

V=  D  a  sinc(A  a/2)exp(-i(k-  — a  +  A  a/2)) 
m  x  m  x 

0) . 

W=  A  a  sinc(A  a/2)exp(i(ic  —a  -  A  a/2)) 
m  x  m  x 

+  “+  + 

Y=  c  a  sine  (A  a/2)exp(i(ic  — a  -  A  a/2)) 
m  x  m  x 

and 

a)  exp(iA  z)-l 

(3-17)  Z  =  exp ( iXjj^z )  [ - - -  ] 

z 

Therefore,  a  typical  term  of  Equation  (3-11), 
integrating  with  respect  to  x'  and  z'  and  fixing  m , 
s ,  is  given  by  : 


(3-18) 


“+  w.  u>?  0)  *  0)  * 

V  <X’2)  *  -  E  Z  «,<=.*♦  +  C,  <=  *  XJ 

*c2 


where  *  indicates  a  complex  conjugate. 
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Chapter  IV 
NUMERICAL  ANALYSIS 


4.1  Determination  of  Axial  Wavenumber 

As  discussed  in  Chapter  II,  eigenvalues  of  the  layered 
waveguide  system  under  consideration  are  determined  by 
examining  the  poles  of  the  Green's  functions,  G^(x|x',z'), 
in  wavenumber  space.  The  real  poles  Xm>  obtained  by  setting 
'H  X  )=0,  are  simply  the  real  roots  of  Equation  (2-26).  As 
stated  in  Section  2.4,  the  axial  wavenumber  must  lie  in  the 
range  k^X^kj  ^or  guided  modes  to  propagate.  Due  to  the 
periodicity  of  the  tangent  function,  several  different 
values  of  the  form  [see  Equation  (2-26)]: 

G  cos  (k-  a)  =Hsin(<„  a) 

^m  Zm  * 

G  =  p  2P  3*  2m*  lm  +  p  1P  2*  Zm*  3m 
H  =  PlP3K:2m  ~p2  ’‘lm’^m 

It  is  then  obvious  from  Equation  (4-1)  that  the  first  real 
root  occurs  for  <  a  in  the  range  0<  ic  a<iT,  the  second  for 


(4-1) 

wh  ere 


-  ^m-"1  *  and  SO  °n*  In  general,  then,  the  axial 
wavenumber  for  guided  mode  m  represented  by  the  positive 
real  roots  of  Equation  (4-1),  Xm>  must  lie  in  the  range 


V-[  fl2 1  x.  i\/^-ra2 


Z  D 

Since  for  guided  modes,  the  lowest  possible 

frequency  that  may  propagate  in  mode  m  is  given  by: 


(4-3) 


f  =  (n-1)  _i, 

"c  ^a  cT^ 
xn  Z 


where  f^  denotes  cutoff  frequency  of  the  m  mode, 
m 

The  asymptotes  described  by  the  left-hand  side  of 
Equation  (4-2)  represent  the  limiting  cases  of  "rigid"  and 
"soft"  (pressure  release)  boundaries  on  medium  II  since  the 
eigenvalues  of  the  systems  are  defined  as 


mw  „  27 

K2m=  T  m  =  0,1, 2, 3,... 


Symmetric  and  asymmetric  modes  for  each  condition  are 
accounted  for  by  m  being  even  or  odd  (e.g.,  for  rigid  walls. 


A  listing  of  the  FORTRAN  program  written  to 
numerically  determine  the  roots  of  Equation  (4-1)  is  given 
in  Appendix  D. 

26 

See  Sec  tion  2.4. 


The  m=0,  or  plane-wave,  mode  can  occur  only  for  rigid 
boundaries  (see  Appendix  B). 


5 


m  even  implies  symmetric  modes;  for  soft  walls,  m  even 
represents  asymmetric  modes).  Figures  7  and  8  show  the 
variation  of  the  axial  wavenumber,  Xm>  with  frequency  for 
the  first  few  modes  of  two,  three-layer  waveguides  with 
semi-hard  (Figure  7)  and  semi-soft  (Figure  8)  boundaries. 
The  asymptotes  given  by  Equation  (4-2)  are  shown  to 
represent  perfectly  rigid  or  soft  walls.  Since  it  was 
assumed  that  c^c^c^  >  the  cutoff  frequency  [Equation  (4-3)1 
is  independent  of  the  speed  of  sound  in  medium  III. 

As  shown  in  Figures  7  and  8,  low  frequencies  tend  to 
propagate  at  or  near  the  phase  velocity  of  the  outer  medium. 
With  increasing  frequency,  the  phase  velocity  asymptotically 
approaches  the  bulk  speed  of  sound  in  medium  II.  This  is 
the  result  of  what  is  referred  to  in  Section  2.8  as 
boundary-induced  dispersion.  The  dispersivity  exhibited  by 
an  acoustic,  slow  waveguide  dramatically  influences  the 
behavior  of  nonlinearly  generated  spectral  components.  A 
discussion  of  these  effects  is  presented  following  the  next 
section . 


4.2  Mod  e  Sha  pe  s 


As  with 

any  bounded 

s  y s  t  em , 

each  mode 

has 

its 

characteristic 

mode  shape 

which 

is  defined 

by 

the 

eigenfunctions  of  the  system.  The  eigenmode  expansions  of 
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Figure  8.  Dispersion  relationship  for  first  few  guided  modes  in  a 

0.1  m. -thick  waveguide  bounded  by  semi-soft  media.  Dashed 
lines  indicate  cutoff  frequencies  for  modes  m=2,3,  and  4, 
and  curves  marked  "asymptotes"  are  bounds  for  each 
mode  (see  Equation  (4-2)). 
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the  velocity  potential  given  in  Chapter  III,  via  Equations 
(3-4,5),  may  be  re-expressed  in  terms  of  cosine  and  sine 
functions: 

CO 

(4-5)  $  =  E  [A^cos^^x)  +  BU1sin(K^x)  ]exp(ix“1z) 

v  a),  q  2q  c  iq  q 

1  q=-°» 

and 

(4-6)  <j>  =  E  [Ca>2cos(ic“2  )  +D<^sin(Klf|c)]exp(ix“22)  . 

0).  S  45  S  L  S  5 

2  s=-°° 

The  velocity  or  pressure  distribution  of  the  source  at  the 
face  of  the  waveguide  (i.e.,  at  z  =  0)  determines  the  extent 
to  which  each  mode  is  excited.  If,  for  example,  the 
velocity  profile  of  a  transducer  matched  that  of  any  single 
mode,  then  only  that  mode  would  propagate  down  the  guide. 
The  series  expansion  of  Equation  (4-5)  would  then  be 
expressed  as  a  single  term  representing  the  particular 
excited  mode,  while  the  coefficients  A and  for  all 

q  q 

other  modes  would  be  identically  zero.  If,  however,  the  end 

conditions  are  not  perfectly  matched,  then  several  modes  may 

be  excited,  all  of  whose  cutoff  frequencies  are  below  the 

2  8 

excitation  frequency. 

Since  the  transverse  pressure  distribution  in  medium  II 
at  a  fixed  range  is  proportional  to  the  time  derivative  of 


2  8 

M.  Redwood,  Mechanical  Wavegul d  e  s  (New  York: 
Pergatnon  Press,  1960),  pp.  77-84  . 
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the  velocity  potential,  then  the  variation  with  x,  or  the 
mode  shape,  is  given  by  the  transverse  component  of  Equation 
(4-5)  for  a  particular  mode  q=m: 

(4-7)  F  (x)  =  A  cos(k  x)  +  B  sin(ic_  x) 

m  m  2m  m  2m  ’ 

where  coefficients  A  and  B  are  determined  by  the  source 

mm 

distribution.  However,  the  ratio  B  /A  may  be  found  through 

mm  ' 

application  of  boundary  conditions  to  Equation  (4-5)  for 
guided  modes : 


(4-8) 


Thus,  the 


B 

m 

A 

m 

function 


P2*3m 

P3K2m 

Fm(x)’ 


normalized  to  A  ,  becomes: 
m 


(4-9) 


Vx)  .  ,  .  Vs™ 

~r~  '  cos(k2»j<)  + 

m  3  2m 


Si„<K2m*) 


The  coefficient  A  outside  the  brackets  has  no  effect  on  the 

m 

shape  of  individual  modes  and  thus  may  be  ignored. 

The  shape  of  the  first  three  modes  at  various 
frequencies  for  a  two-layered  medium  is  depicted  in  Figures 
9-11.  These  cases  correspond  to  a  large  flat  plate  of 
silicone  rubber  immersed  in  water. 
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Figure  9.  Transverse  pressure  distribution  in  medium  II  at  fixed 
range  for  m=l  mode  at  5,  10,  and  15  kHz.  Vertical  axis 
represents  vertical  distance  from  medium  III  and  horizontal 
axis  is  given  by  Equation  (4-9).  Physical  parameters 
correspond  to  a  0.1m. -thick  slab  of  silicone  rubber  immersed 
in  water.  SI  units  used. 
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Figure  10.  Transverse  pressure  distribution  in  medium  II  at  fixed 

range  for  m=2  mode  at  10,  15,  and  20  kHz.  Vertical  axis 
represents  vertical  distance  from  medium  III  and  horizontal 
axis  is  given  by  Equation  (4-9) .  Physical  parameters 
correspond  to  a  0.1m. -thick  slab  of  silicone  rubber  immersed 
in  water.  SI  units  used. 


MODE  m=3 
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Figure  11.  Transverse  pressure  distribution  in  medium  II  at  fixed 

range  for  m=3  mode  at  15,  20,  and  25  kHz.  Vertical  axis 
represents  vertical  distance  from  medium  III  and  horizontal 
axis  is  given  by  Equation  (4-9).  Physical  parameters 
correspond  to  a  0.1m. -thick  slab  of  silicone  rubber  immersed 
in  water.  SI  units  used. 
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Hence,  the  phase  velocity  and  density  correspond  to  previous 

29 

experimental  work.  It  should  be  noted  that  the  increase  in 
amplitude  at  higher  frequencies  for  each  mode  further 
demonstrates  the  effect  of  dispersion  discussed  above. 

Finally,  the  variation  of  mode  shape  with  changing 
boundary  conditions  for  the  lowest  mode  (m«l)  is  shown  in 
Figure  12.  The  appearance  of  spreading  seen  for  lower 
values  of  pc  (i.e.,  softer  boundary  conditions)  at  a  given 
frequency  is  explained  by  Figures  7  and  8.  The  harder  the 
boundaries,  the  less  influence  media  I  and  III  (which  in  all 
of  the  cases  depicted  are  identical)  have  on  the  acoustic 
field  in  medium  II. 


4.3  Behavior  of  Nonlinearly  Generated  Components 
in  Dispersive  Media 

The  lowest  spectral  component  produced  by  the  nonlinear 
interaction  of  two  primary  waves  of  frequencies  a)  ^  and  is 
the  difference  frequency  u)  (=u^-u)2)*  The  propagation 
constant  of  the  volume  source  produced  via  interaction  of 
the  primary  waves  at  this  frequency  is  determined  by  a 


29 

J.D.  Ryder,  P.H.  Rogers  and  J.  Jarzynski,  "Radiation 
of  Difference-Frequency  Sound  Generated  by  Nonlinear 
Interaction  in  a  Silicone  Rubber  Cylinder,"  Journal  of  the 
Acou.stlcal  Society  o  f  Amer  ic  a .  59  (  1  976),  1077-1086;  See 
»  so  Figure  6. 
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(a)  (b)  (c) 


MEDIUM 

PHASE  VEL. 

DENSITY 

PHASE  VEL. 

DENSITY 

PHASE  VEL. 

DENSITY 

I 

II 

III 

1500.0 

1000.0 

1500.0 

1000.0 

200.0 

1000.0 

1500.0 

1000.0 

1500.0 

1000.0 

1000.0 

1000.0 

1500.0 

1000.0 

1500.0 

1000.0 

5000.0 

1000.0 

MODE  m=l  THICKNESS  =  6.1  m.  SI  UNITS 

FREQ. (kHz)  =8.0 

CO 


FjCXD/A  j 


Figure  12.  Transverse  pressure  distribution  in  medium  II  for  m=l  mode 
at  8.0  kHz  for  a  0.1m. -thick  slab  of  varying  density 
immersed  in  water.  Axes  are  same  as  for  Figure  9. 
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vector  combination  of  the  primary  wavenumbers*  If  this 
vector  combination  does  not  correspond  to  the  natural 
propagation  constant  of  the  difference  frequency  in  the 
medium,  then  asynchronous  interaction  will  occur,  resulting 
in  the  "spatial  beating"  effect  depicted  in  Figure  13  for 
Az>0.  If, however,  the  primary  waves  which  give  rise  to  the 
difference-frequency  component  propagate  at  the  same  phase 
velocity  as  the  latter,  then  Figure  13  shows  that 
synchronous  interaction  or  "spatial  resonance"  occurs 
(Figure  13,  Az=0).  Under  this  condition,  the  amplitude  of 

the  nonlinearly  generated  component  will  grow  monotonically 
with  range  until  enough  energy  has  been  transferred  from  the 
primary  waves  to  significantly  reduce  the  strength  of  the 
volume  source  distribution  (i.e.,  forcing  function). 

In  general,  the  vector  A  is  used  to  represent  the 
wavenumber  difference  between  the  primary  waves  and 
nonlinearly  generated  frequency  components.  Thus,  for  the 
case  of  the  sum-  or  difference-frequency  components 
resulting  from  a  bifrequency  primary  wave  interaction, 

(4-10)  A  =  nk.±  rak.  -  k 

1  i  nm 

In  an  unbounded,  dispersionless  medium,  the  resonance 
condition  can  be  satisfied  only  if  the  wavefronts  of  the  two 


waves  are  parallel  and  both  propagate  in  the  same  direction. 
Under  these  conditions,  the  linear  relationship  between 
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(9 


RANGE  CM} 


Figure  13.  Magnitude  of  axial  component  of  velocity  potential,  Z,  in 
medium  II  as  given  by  Equation  (3-17)  for  several  values 
of  real  Az  (=1,2, 4, 8),  hence  damping  is  neglected  (i.e.  a=0) . 


wavenumber  and  frequency  (k-^/c)  necessitates  synchronous 


interaction*  However,  in  a  dispersive  medium  such  as  a 
waveguide,  resonant  interaction  is  not  so  straightforward 
and  will  be  treated  in  Section  4.4* 

As  previously  stated,  the  analytical  basis  for  Figure 
13  may  be  found  by  examining  the  axial  component  of  the 
velocity  potential  for  the  difference  frequency  in  medium 
II,  which  was  derived  in  Chapter  III  [i.e.,  Equation 
(3-17)].  Moreover,  since  acoustic  waves  in  medium  II  may  be 
characterized  by  their  axial  wavenumber  for  each  mode,  the 
vector  A  becomes  a  scalar  quantity  defined  for  the  sum  or 
difference  frequency  by 


(4-11) 


A  =  x^l  +  x“2  ~ 
z  Aq  —  s 


where  the  integer  subscripts  q,  s,  and  m  refer  to  the  modes 
in  which  the  primary  frequencies  (  and  )  and  the  sum  or 
difference  frequencies  are  respectively  under  consideration. 
The  variation  of  the  velocity  potential  with  range  is 
therefore  described  by  the  following  equation  from  Chapter 
III: 


uj+  exp(iA  z)-l 

(3-17)  Z  =  exp(ixm-z)[ - -  1 

z 


jn urns'. 
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where  the  magnitude  of  Z  is  plotted  in  Figure  13  for  several 

real  values  of  A  •  Notice  that  as  A  Increases,  the  maximum 
z  z 

allowable  amplitude  of  the  sum-  or  difference-frequency 

velocity  potential  decreases.  In  fact,  the  maximum 

allowable  amplitudes  |Z|  occur  at  distances  z  from  the 

max  m 

source,  where 


(4-12) 


lZLax  *  ^  z=z  A 
m  z 


(4-13)  2=7“  n-1.3,5 . 

m  A 

z 

Thus,  in  order  to  transfer  as  much  energy  as  possible  from 

the  primaries  into  a  nonlinearly  generated  frequency 

component,  such  as  the  sum  or  difference  frequency,  A^  must 

be  minimized  for  that  frequency. 

It  can  be  seen  from  Equation  (4-12)  that  the  maximum 

amplitude  of  the  nonlinearly  generated  component  can  become 

extremely  large  for  small  A  •  In  fact,  for  A  =0 

z  z 

(resonance),  the  amplitude  grows  linearly  with  range  as 
shown  by  Figure  13.  This  apparently  unrealistic  situation 
may  be  resolved  by  considering  the  fact  that  the  growth  of  a 
nonlinearly  generated  component  is  accompanied  by  the  decay 
of  the  fundamental  waves.  Eventually,  the  decreased 
strength  of  the  primary  fields  no  longer  permits  a  nonlinear 


interac  tion 


•  A-. 
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All  real  physical  systems  involve  a  certain  amount  of 

damping,  which  is  accounted  for  by  including  an  imaginary 

30 

term  in  the  axial  wavenumbers.  In  this  instance.  Equation 
(3-17)  becomes: 


(4-14) 


1  -exp(iA  z)  exp(-CLz) 

Z  =  exp(iX  ^z)  exp  (-a  z)  - rr - 

r  m  r  la,  am  -  iA 


where  a  and  a  represent  damping.  The  magnitude  of  Z 

T  0)±  * 

given  by  Equation  (4-14)  versus  range  for  fixed  Az>  with 
absorption  as  the  parameter,  is  illustrated  in  Figure  14. 


4.4  Dispersion  in  Medium  II 

As  established  in  the  previous  section,  the  maximum 

amount  of  energy  transferred  to  any  nonlinearly  generated 

frequency  component  in  a  lossless,  dispersive  medium  is 

3  1 

determined  by  the  difference  in  axial  wavenumbers  A  .  The 

term  A  as  described  by  Equation  (4-11)  is  a  function  of 
z 

several  independent  variables.  The  \  versus 

m 

relationships  (e.g.,  Figures  7  and  8)  which  characterize  the 
dispersivity  of  the  waveguide  are  determined  by  the  physical 


30 

Valid  only  for  "weak"  nonlinear  interactions. 
^See  Equations  (4-11)  and  (4-12). 


RANGE  CM} 


Figure  14 


Magnitude  of  .axial  component  of  velocity  potential, 
in  medium  II  as  given  by  Equation  (A-1A)  for  Az=4.0 
increasing  amounts  of  absorption  (i.e.,  a=.001,.05, 
SI  units  used. 


Z, 

with 

.5,1.0) 
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parameters  of  the  media  (i.e.,  density  and  bulk  phase 
velocity  of  media  I,  II,  and  III,  and  thickness  of  medium 
II).  Once  this  relationship  is  defined,  &z  may  be 
determined  for  the  sum-  or  difference-frequency  component  of 
any  given  pair  of  primary  waves,  as  long  as  it  is  known  in 
which  modes  all  four  frequencies  propagate. 


The  variation 

o  f 

A  for  the 

z 

nonl inea  rly 

generated 

difference 

frequency 

u> 

for  a  fixed 

ratio  W1/(D_  is 

shown  in 

Figure  15. 

Again , 

the 

physical  parameters  are 

chosen  to 

correspond 

with  previous  empirical 

data  (i.e., 

silastic 

rubber  in 

seawater) 

32 

• 

The  primary 

waves  as  well  as  the 

difference 

frequency 

are 

assumed  to 

propagate  in 

the  first 

mo  de  . 


32see  e.g.,  Ryder,  Rogers,  and  Jarzynski,  pp. 
1077-1086. 


MEDIUM  PHASE  VEL.  DENSITY  THICKNESS 


PRIMARY  MODE  q,s  =  l 


I 

II 

III 


1500.0 

1000.0 

1500.0 


1000.0 

1000.0 

1000.0 


00 

varies 


DIFF.  FREQ.  MODE  m=l 
FREQ.  RATIO  f  /f_ =  10 
THICKNESS*  =  .05,  .1,  .2 


SI  UNITS 


O 


Figure  15.  Variation  of  h?  with  difference  frequency  (f_)  for  fixed 
ratio  of  primary  frequency  (fj)  to  f_  (i.e.,  f^/f_=10.0). 
Phase  velocities  and  densities  correspond  to  silicone  rubber 
immersed  in  water.  It  is  assumed  that  the  difference 
frequency  as  well  as  the  primaries  propagate  in  mode  m=l. 

The  three  curves  are  for  thicknesses  of  .05,  .10,  and 
.20  meters. 


Chapter  V 


CONCLUSIONS 


A  theoretical  investigation  of  weak,  nonlinear  acoustic 

wave  interactions  in  a  three-layered  medium  has  been 

presented  in  this  thesis.  Under  the  conditions  explained  in 

Chapter  I,  an  acoustic,  slow  waveguide  has  been  analyzed 

both  numerically  and  theoretically,  and  various 

characteristics  of  linear  waveguide  theory,  as  well  as 

finite-amplitude  phenomena,  have  been  treated. 

The  dispersivity  for  guided  mode  propagation  in  an 

acoustic,  slow  waveguide  is  clearly  defined  by  the 

characteristic  equation  derived  in  Chapter  II.  With  the  aid 

of  the  computer  program  listed  in  Appendix  D,  exact 

numerical  data  can  be  calculated  given  the  physical 

parameters  of  the  system.  Once  the  dispersion  relationship 

is  established,  the  maximum  allowable  amplitude  for  any 

nonlinearly  generated  component  may  be  calculated  simply  by 

determining  the  corresponding  value  of  A  . 

z 

Expressions  which  may  be  used  to  determine  the 
••'inversion  efficiency  of  weak  parametric  interactions  in 
.mple,  layered  media  (i.e.,  liquid-like  media)  have  thus 
i>vrloped.  Moreover,  a  comparison  of  conversion- 
•  ■i  v  pnhancement  in  a  waveguide,  relative  to  that  in 


^  IS 
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an  unbounded  medium,  can  now  be  realized.  However,  in  order 
to  place  an  upper  bound  on  conversion-efficiency  enhancement 
via  boundary-induced  dispersion,  a  complete  numerical 
analysis  of  strong  wave  interactions  is  required.  Such  an 
analysis,  which  would  involve  implementation  via  digital 
computer  solution  of  many  coupled  partial  differential 
equations,  is  outside  the  scope  of  this  investigation.  It 
is  therefore  recommended  for  future  study. 


APPENDIX  A 


Determination  of 


3V(x) 


3x 


X=X. 


m 


From  Equation  (2-21): 
(A-l) 


where 


'f(x)  =  (a+d)sin(><2a)  +  i(b+c)cos(<2a) 


a  =  p2  k'1K3 


b  =  plp2  K2K3 


c  =  p?  p3k3k2 


d  P1P2P3IC2 


Then, the  derivative  is  given  by: 

(A-2)  T’  (X  ) 

m 


3T 

3X 

X=Xm 

3k . 

_ ^  _1_ 

^.2  + 

3k-. 

3ki 

3X 

3k2  3x 

3X 

Also  , 

(A-3) 


3X 


=  _  JL 


K  , 
1 


i  =  1,2,3 


and 


3T 

3^ 


P.  k  sin(K  a)  +  ip  p  K„cos(k-  a) 


(A-A) 


d< 


3 


Therefore 
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Appendix  B 


PLANE-WAVE  MODE 


As  discussed  in  Chapter  II,  the  Green's  functions  for 
the  waveguide  are  expressed  as  a  summation  of  the  residues 
of  the  transformed  Green's  functions,  G^ ,  at  the  real  poles 

X  of  G  .  The  real  roots  of  Equation  (2-26),  restated  here, 

X 

represent  these  poles: 


(B-l) 


tan(<2ma) 


<lm'C2mP2Pj  *  IC2m<3mPlP2 

2  2 
IC_  P,P,  -  K,  1C-  p_ 

2m  13  lm  3m  2 


One  root  in  particular  of  this  equation,  v  =k_  (or  k  .  *=0), 

m  2  2m 

is  not  consistent  with  physical  or  mathematical  assumptions. 

It  is  the  purpose  of  this  section  to  examine  the  behavior  of 

G  at  the  root  Xra=^2  anc*  to  determine  its  physical 

X 

sign  if icance  . 

Assume,  as  in  Equation  (2-27),  that  G  is  expressible 

X 

in  the  form 


(B-2)  G^(x  |  x’ ,  z' )  -  exp(ix(z-z')) 

where  H'(y)  is  given  in  Equation  (2-21),  and  F(X),  which 
changes  with  each  region,  may  be  deduced  from  Equations 
( 2- 1 5 ) - ( 2- 1  8 )  .  For  region  IIA  (i.e.,  a>x>x'); 
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(B  3)  F(x)  -  Aexp(-i(0-j>))  ±  Bexp  (i  (SH-j>) )  +  Cexp  (-!(£>+$) )  +  Dexp(i(0-<t>)) 

4x2 

0  -  «2(x-a) 

where 

<?  *  <2*' 


and  coefficients  A,  B,  C,  and  D  are  given  by  Equation 
(2-20).  Then,  since 


(B-4) 


limit 


is  indeterminant,  L'Hopital's  Rule  may  be  applied  with  the 
following  results: 


(B-5) 


limit  G  =  limit 
X 


X 


X  +  K 


3F(x) 

3X 


x  =  xTO  3  K 

m  i 


3y(x) 

3X 


X  =  X  3  k7 
m  z 


-k,  (x-a)+<_  x 1  —  x ,  (x'(x-a))+l 

lm _ 3m _ lm  3m 

k.  +<_  +x.  a 

lm  3ra  lm  3m 


Thus,  for  guided  modes,  and  m  are  both  positive  real 

values;  therefore,  the  transformed  Green's  function  Gy  in 
region  IIA  remains  finite  at  Xm=k2*  Hence,  k2*X  is  not  a 
oingularity  of  G  even  though  it  is  a  root  of  Equation 

A 

(2-26).  A  similar  process  verifies  that  for  all  regions 
(i.e.,  I,  IIA,  IIB,  and  III), 


(B-6)' 


limit  G  (x | x' ) 

<2*  0  L  X 


remains  finite*  Physically,  this  means  that  an  incoming 
wave  propagating  purely  in  the  axial  direction  of  the 
waveguide  (i.e.,  X  implies  a  planar  wavefront)  will  not 
propagate  as  a  trapped  wave,  but  will  decay  as  energy  is 
radiated  to  the  outer  media.  If,  however,  medium  II  had 
rigid  boundaries,  then  the  plane-  wave  mode  would  propagate 
as  a  guided  mode. 
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Appendix  C 

GREEN'S  FUNCTIONS  FOR  MEDIA  I  AND  III 

Expressions  derived  in  Chapter  II  for  G*  represent 
Green's  functions  for  medium  II  only.  That  is,  it  was 
assumed  that  the  point  source  was  in  the  waveguide  (recall 
that  0<x'<a).  If,  however,  sources  occurred  outside  the 
guide,  or  energy  originally  radiated  from  the  guide 


penetrated  back  into  medium  II  via  some  coupling  mechanism 
(e.g.,  corregated  surface  structure),  then  Green's  functions 
for  media  I  and  III  would  be  required  in  order  to  evaluate 
these  effects. 

As  discussed  in  Chapter  II,  the  transformed  Green's 

functions  G  ,  for  the  waveguide,  were  derived  for  each 

X 

region  and  are  given  by  Equations  ( 2- 1 5  ) - ( 2- 1 8 ) .  In  order 
to  distinguish  these  from  the  following  expressions  which 
are  obtained  by  assuming  the  point  source  is  in  medium  I  or 
III,  a  superscript  will  be  used  to  denote  medium  number 

(i.e.,  G1  and  G *** correspond  to  the  point  source  in  media  I 

X  X 

and  III,  respectively). 


The  matrix  equations  whose  solutions  determine  the 
transformed  Green's  functions  for  media  I  and  III  are  given 


by  Equations  (C-l)  and  (C-2),  respectively. 

Following  are  expressions  for  the  transformed  Green's 
functions  for  medium  I  (i.e.,  x'>a): 


Region  IA  ,  x>x ' : 


(C-3  )  G  1  =  exp(-ixz’ iexpCiK^)  [ 


exp (i<2a) [Cexp (iK^(x' -a) )  + 

Hx) 


Dexp(-itc1(^-a))]+  exp(-iK9a)  [Eexp(i<1(x,-a))  +  Fexp(-i<1(x'-a))  ] 


'Hx) 


Region  IB,  x'>x>a: 

exp  ( i* 9 a )  [  Cexp  ( iic .  (x-a ) )  + 

(C-4)  G^  *  expC-ixz^expCiK^')  [ - - - - 


Dexp (-iK 1 (x-a) ) ]  +  exp (-ii^a) [ Eexp ( i<1 (x-a) ) +  Fexp (-ix^  (x-a) ) ] 

~~?0c) 


Reg  ion  II,  a  >x>0 : 


t  P1P^<1ic2cos(k2x)  -  ip^p2<^ic3sin(’c2x^ 

(C-5 )  g  1  =  4exp(-ixz')exp(iic1x')  - “ 


Reg  ion  III,  x<0: 


G 

X 


I 


P1P9<1<?4exp(-ixz'  )exp(-iK;3x)exp(iK1x') 

nx) 


(C-6) 


8 

The  transformed  Green's  functions  for  medium  III  (i.e,  x'<0 
are : 


Region  I,  x>a: 


(C-7) 


III  4p 2^ 3lc2K:3e3<P ^-iXZ *  )exp(-iK3x,)exp(iK]x) 

Gx  =  nx) 


Region  II,  a>x>0: 


(C-8) 


4exp(-ixz '  )exp(i<1  a)exp  (-ii^x1 )  P2P 3IC^<3C0S (x-a) )  + 

[  nx) 


iplP3K2K3Sin^2^x_a^ 

?(x) 


Region  IIIA,  x'<x<0: 


(C - 9 ) 


III 


=  exp(-ixz'  ^xpUK^expC-i^x’)  [ 


exp(i<3x) [Dexp(ix?a)  + 

nx) 


Fexp(-iK2a)  ]  +  exp (-i^x)  [Gexp (ii<2a)  +  HexpC-iic^a)  ] 
—  — 


1 


Region  1 1 1  B  ,  x  <x  '  : 


.  TTT  exp(k.a)[Dexp(k.x')  + 

(C-l  0)  G^iil  =  eXp(_iX2-)exp(ix1a)exp(-U-3x)  [ - - =? - 


Gexp(-iK_x' )  ]  +  exp(-i«.a)  (Fexp(i»c-x'  )  +  Hexp (-ix.x' )  ] 


The  Green's  functions  are  then  obtained  via  an  inverse 
transformation  involving  complex  integration  techniques 
(i.e.,  residue  theory)  as  outlined  in  Chapter  II.  Following 
are  the  results  for  media  I  and  III. 

Green's  functions  for  the  mth  forward-guided  mode  for 
medium  I  are: 


Region  IA,  x>x' 


t  exp(i<  a) [C  exp(-<  (x’-a))  + 

(C-H)  Gm  =27Tiexp(iXra(z-Z’))exp(-Klnx)  [ -  ) - “ - 


PmeXp(Klm(x,~a))J  +  exp(-i<2nia)[Emexp(-lclm(x’-a)) 

V'(X  ) 

m 


Fmexp(Klm(x,~a))] 

m 


Region  IB,  x'>x>a: 


(c  I  exp(x  (a-x))[C<Ttexp(ix  a)  + 

(  Gm  =  2Triexp(ixm(z-Z’))exp(-<lmx')[ - i - - =SL_ 

m 


Emexp(-i<2ma)]+exp(-xlm(a-x))[DmeXp(ix2ma)  +  Fmexp  ("i,C2ma)  ] . 


v’  (xm) 


«SHS«|.>..W  <v»  V.,  M HA.-*..*  ■■ ■ ■“ 
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Region  II,  a>x>0: 


(C-13) 


P  P  *C  K  COS(K  X)  + 

Gm  =  -Sirexp  ( ’ >  >  exp  (-k^x  ’ )  [ - } - 


m 


V’(X  ) 
m 


Region  III,  x<0: 


(C-14) 


G  = 
a 


"8irp  lp  (iXm(z“Z  }  eXp (K3mX>  exp  (“XlmX ' } 

V’  (X  ) 
m 


And  the  Green's  functions  for  medium  III  are: 


Region  I,  x>a: 


(C-15) 


III  ~8"P2P3K2m<3mexp(ixm(2"2,))exp(<3mX,)exp(~<lmx) 


m 


r  (X  ) 

m 


Region  II,  a>x>0: 


Region  IIIA,  x'<x<0: 
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(C-17)  G  =  27riexp(iX  (z-z'))exp(-<  a)exp(<-  x')  x 

m  m  Im  jm 


exp(-K,  x) (D  exp(i<-  a)  +  F  exp(-i<_  a)] 
_ 3m _ m  _ Zm _ m  Zm 

r(X  ) 

ra 

exp«3nx)  [Gmexp(iK2ma)  +Hmexp(-i<2ma)] 

r  (X  ) 

m 


Region  IIIB,  x<x': 

(C-18)  G  111  =  2iriexp(iX  (z-z ’ ) ) exp (-£.  a)exp(ic  x)  x 

m  m  iiti  jm 

exp  ( i<:ma )  [  Dmexp  (-<3mx ' )  +  G^exp  (k^)] 

[  rex  ) 

m 

exp(-itc_  a)[F  exp(-ic_  x’)+H  exp(ic-  x1)] 

_ £ra _ ITT _ Jm _ m _ Jm  , 

v  (x  )  ] 

m 

where  the  following  substitutions  have  been  made  in 
Equations  (C-3)-  (C-10); 


(C-19)  f  (x)  =  2iexp(iic  a)  [Aexp(iV:2a)  +  Bexp(-iK'2a)  ) 


(C-20) 


a  p0p3<^x2 
2 

b  =*  p  p  *2 
c  -  P22<1x3 
d  *  P  ^P  9  ^  'y  ^ 
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(C-21)  A— a+b+c-d 


B  =  -a -b -c  -  d 


C-a+b-c-d 


D=a-b-c+d 


E-a-b+c-d 

F-a+b+c+d 

G--a+b-c+d 

H»-a-b+c+d 


Equations  ( C  —  1 1 )  through  (C-18)  represent,  as  in  Chapter  II, 
the  mt^1  guided  mode  in  an  eigenmode  expansion  of  the  Green's 
functions*  The  poles  of  are  found  by  setting  Equation 
(C-19),  ^ (x) —0 •  Then,  for  guided  modes 


(C-22) 


K1  “  1Klm 

k2  =  K2m 
k3  *  1K3m 


and  the  terms  given  by  Equation  (C-20)  become: 

<C-23> 

2 


ip2P3KlmK2m 


m 


PlP3K2m 


Cm  "  _P2  KlmK3m 
dm  *  lplp2K2mK3m  • 

The  coefficients  A  ,B  ,...H  are  the  same  functions  as  given 

mm  m 

in  Equation  (C-21)  except  for  the  use  of  the  subscripted 
variables  defined  in  Equation  (C-23)  in  place  of  the  a,b,c,d 
found  in  Equation  (C-20)* 

Notice  that  the  determinants  of  the  matrices  of 
Equations  ( C  —  1  )  and  ( C  —  2 )  are  equal  [Equation  (C-19)),  but 


" '  * woinwiS'Vw'.T  ■ 
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differ  from  the  determinant  found  in  Chapter  II  [Equation 
(2-21)].  The  derivative  of  Equation  (C-19)  evaluated  at  the 
poles  X-Xj,,. 


(C-24) 


*'<  Xj 

m 


3X  X  „X 


> 


therefore  should  not  be  confused  with  the  expression  derived 
in  Appendix  A,  Equation  (A-7),  which  was  restated  in  Chapter 
II.  This  expression  [Equation  (2-38)J  is  derived  by 


assuming  the  point  source  lies  in  medium  II. 


Appendix  D 
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NUMERICAL  EVALUATION  OF  EIGENVALUES: 
COMPUTER  PROGRAM  LISTING 


The  FORTRAN  program  named  AXIAL  (Figure  16)  was  written 
by  the  author  for  the  purpose  of  numerically  investigating 
various  characteristics  of  an  acoustic,  slow  waveguide.  The 
eigenvalues  of  the  guided  modes  in  medium  II  are  obtained  by 
evaluating  the  real  roots  of  the  characteristic  equation 


(D-l) 


tan(<2ma) 


P2P3IClm<2m'f~PlP2IC2m'C3m  =  G 
PlP3K2m  "P2  Klm'C3m 


•  t 


In  order  to  avoid  numerical  problems  at  the  discontinuities 

of  the  tangent  function,  a  similar  form  of  Equation  (D-l) 

33 

was  used  in  the  program: 

(D-2)  F(x)  “  Gcos(<2ma)  -  Hsin(K2na) 

AXIAL  utilizes  an  IMSL  subroutine  (ZBRENT)  to  locate 
the  real  roots  of  an  external  function  F(x).  These  real 
roots  are  stored  in  a  multidimensional  array  TRAP(m,n)  where 


33See  also  Equation  (4-1), 
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"m"  denotes  mode  and  "n"  frequency  number*  Since  ZBRENT  was 

designed  to  locate  the  root  x'  of  a  function  between  any  two 

endpoints  a  and  b  such  that  the  product  of  F(a)  and  F(b)  is 

less  than  zero,  care  must  be  taken  to  insure  that  exactly 

one  root  exists  in  that  interval  before  executing  ZBRENT. 

It  follows  from  the  discussion  in  Section  4.1  that  the  n1*1 

real  root  lies  between  the  asymptotes  described  in 

equation  (4-2).  Hence,  the  interval  containing  one  and  only 

one  root  of  F(x)  is  well  defined. 

Once  the  roots  have  been  located  and  stored,  the  array 

TRAP  is  plotted  via  the  CCS  graphics  system  supported  by  The 

34 

Pennsylvania  State  University  Computation  Center. 
Examples  of  the  graphics  output  may  be  seen  in  Figures  6,  7, 

and  8 . 


3  4 

CCS  is  an  abbreviation  for  California  Computer 
Products,  Inc.  (CalComp)  who  originally  developed  the 
sof  tware • 
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AUTHOR: 

DAVID  M.  YEAGER 
ACOUSTICS  DEPARTMENT 
PENNSYLVANIA  STATE  UNIVERSITY 

DATE  WRITTEN: 

AUGUST  22,1979 

PROGRAM  NAME: 

AXIAL 

PURPOSE: 

THE  PURPOSE  OF  THIS  PROGRAM  IS  TO  COMPUTE  THE  POSITIVE,  REAL  ROOTS 
X  OF  AN  EXTERNAL  FUNCTION  F(X)  WHERE  'F'  IN  THIS  CASE  IS  THE 
CHARACTERISTIC  EQUATION  FOR  A  TWO-DIMENSIONAL  FLAT-PLATE  'SLOW' 
WAVEGUIDE.  HENCE  THE  PHASE  VELOCITY  IN  MEDIUM  TWO  IS  LESS  THAN 
IN  MEDIA  1  OR  3.  IN  FACT  IT  IS  ASSUMED  THAT  C3  >  Cl  >  C2. 

THESE  ROOTS  REPRESENT  THE  EIGENVALUES  OF  THE  GUIDE  AND  ARE 
GROUPED  ACCORDING  TO  MODES  BEFORE  BEING  PLOTTED.  THE  VERTICAL 
AXIS  IS  FREQUENCY  AND  WAVENUMBER  IN  THE  AXIAL  DIRECTION  IS 
DISPLAYED  ON  THE  HORIZONTAL  AXIS. 


REAL*8  EPS,A,B,F 
REAL  MAXFR 
CHARACTER*8  DA 

DIMENSION  GAM1 (202) ,GAM2(202) ,GAM3(202) ,FRQ(202) ,TRAP(10,202) , 
CFC(30) .ALPH1 (202) ,ALPH2(202).ALPH3(202) ,DELTA(202) ,C(3) ,R(3) 
COMMON /AREA 1 /GAM II , GAM2I ,GAM3I , RO 1 , R02 , R03, I , THICK 


SEE  LISTING  OF  ZBRENT  FOR  EXPLANATION  OF  'EPS.NSIG.MAXFN' 

'C1.C2 ,C3'  REPRESENT  BULK  PHASE  VELOCITIES  FOR  MEDIA  1,2,3 
RESPECTIVELY.  'THICK'  IS  TOTAL  THICKNESS  OF  MEDIUM  2  SINCE  THE 
BOUNDARIES  ARE  AT  X=0  AND  X-A.  'ROl ,R02,R03'  REPRESENT  DENSITIES 
OF  THE  THREE  MEDIA.  'MAXFR'  IS  THE  MAXIMUM  FREQUENCY  PLOTTED. 


READ (5 , 8000)  EPS.NSIG.MAXFN 

READ (5 ,8001 j  Cl .C2.C3, THICK, MAXFR 

READ (5, 8002)  R0l,R02,R03 

MAXFN 1-MAXFN 

TWOPI“2.*3. 14159 

DUMMY 4 -MAXFR*. 001 

WRITE(6 , 9000)  Cl. C2.C3, THICK 

WRITE (6, 9005)  R0l,R02,R03 

WRITE(6 , 9001 )  DUMMY 4 

WRITE (6 ,9002)  EPS.NSIG.MAXFN 


CALCULATE  CUTOFF  FREQUENCY  FOR  MODE  M,  STORE  IN  FC(M) , 
AND  WRITE  RESULTS  FOR  EACH  MODE  TO  BE  PLOTTED. 


M-2 

FC( 1 )“0 . 

ANUM=2.*(M-1 .) 

DENOM-4 . *THICK*SQRT ( 1 . /C2**2-l . /Cl**2) 


Figure  16.  FORTRAN  listing  of  root-finding  computer  program  AXIAL 
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7000 

C 

C 

C 

C 

30 


C 

C 

C 


c 

c 

c 

c 


50 

C 

C 

C 

40 


C 

C 

c 

c 

10 

c 

c 

c 

c 

c 

c 


FC (M) -ANUM/DENOM 
FCM-FC(M) 

IF  (FCM  .GT.  MAXFR)  GO  TO  30 

WRITE (6, 9004)  M,FCM 

M-M+l 

IF(M  .GT.  10)  GO  TO  7000 
GO  TO  29 
WRITE(6,9020) 

'FREQ'  IS  THE  INITIAL  FREQUENCY,  'P'  IS  THE  FREQUENCY  INCREMENT, 
WRITE  TOTAL  NUMBER  OF  MODES  TO  BE  PLOTTED 

FREQ-100.  ,  _ 

P-(MAXFR-FREQ)/200. 

X-M-l .  ,  v 

NMODE-IFIX(X) 

WRITE (6, 9009)  NMODE 

MAIN  LOOP  ITERATED  FOR  200  FREQUENCIES. 

DO  5  1-1,200 
J-0 

FRQ(I)-FREQ 

GAMl (I)-TW0PI*FREQ/C1 

GAM2 (I)-TW0PI*FREQ/C2 

GAM3(I)-TWOPI*FREQ/C3 

GAMII-CAMI(I) 

GAM2I-GAM2(I) 

GAM3I-GAM3U) 

BELOW  CUTOFF  FREQUENCY  FC{M)  THE  MTH  MODE  IS  INITIALIZED  TO  THE 
WAVENUMBER  IN  MEDIUM  1. . .GAM1I«TW0PI*FREQ/C1 . 

DO  50  N-l, NMODE 
TRAP(N,I)-CAM1I 
CONTINUE 

SELECT  ENDPOINTS  OF  INTERVAL  TO  BE  SENT  TO  ZBRENT. 


J-J+l 

A-GAM1I+.0001 

B-GAM2I-.0001 


IF(J  .GT.  1)  B-SQRT(GAM2I**2-(TWOPI*(J-1.)*2./(4.*THICK))**2)-.01 
IF  (FREQ  .GT.  FC(J+1))  A«SQRT(GAM2I**2-(TWOPI*(J*2.  )/(4.*THICK) 

\  \  i  r\  f 


C)**2)+.01 
FA-F(A) 

IF(FA*^B  .LT.  0.)  GO  TO  10 


IF  F(A)  AND  F(B)  HAVE  THE  SAME  SIGN  THEN  NO  ROOT  EXISTS  IN  THAT 
INTERVAL.  GO  ON  TO  NEXT  MODE  OR  FREQUENCY. 


GO  TO  41 
MAXFN-MAXFN1 


'ZBRENT'  IS  AN  IMSL  SUBROUTINE  WHICH  LOCATES  THE  ROOT  OF  AN 
EXTERNAL  FUNCTION  'F(X)'  BETWEEN  ENDPOINTS  'A'  AND  'B'. 

THE  ROOT  IS  RETURNED  AS  THE  VALUE  'B'  SUCH  THAT  'F(B)-0'. 


Figure  16.  (continued) 
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C 

C 

C 

41 

20 

4 

5 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 


C 

C 

C 


C 

C 


C 

C 


60 

C 

55 

C 

C 

C 


CALL  ZBRENT(F,EPS,NSIG,A.B,MAXFN,IER) 

IF(IER  .EQ.  129)  CO  TO  20 

STORE  ROOTS  FOR  MODE  J,  FREQUENCY  I  IN  'TRAP'. 

TRAP(J,I)-B 

IF(FREQ  .GT.  FC(J+1))  GO  TO  40 
GO  TO  4 

WRITE (6, 9006)  MAXFN.I 

FREQ-FREQ+P 

CONTINUE 


BEGIN  PLOTTING  FREQUENCY  VS.  AXIAL  WAVENUMBER  VIA  'CALCOMP' 
GRAPHICS  SUBROUTINES . 

PLOT  GAM2  FIRST  SINCE  IT  HAS  THE  LARGEST  WINDOW.  THEN  PLOT 
GAM1  AND  GAM3.  USE  SCALE  AND  TRANSLATION  FACTORS  STORED  IN 
GAM2(202)  AND  GAM2(201). 


CALL  PLTTYP (4662 ,6,7) 

CALL  START 

MOVE  ORIGIN  4.0  IN.  OVER  AND  2.0  IN.  UP  FROM  LOWER  LEFT  CORNER. 

CALL  PLOT(4. 0,2. 0,-3) 

CALL  SCALE(GAM2,9. 0,200.1) 

CALL  SCALE (FRQ. o . 5, 200,1) 

CALL  AXIS(0. 0,0.0,  WAVENUMBER 
CGAM2(20l),GAM2(202)) 

CALL  AXIS ( 0 . ,0. , ' FREQUENCY' , 9 , 6 . 5, 90 . , FRQ(201 ) , FRQ (202 ) ) 

CALL  DATE (DA) 

CALL  SYMBOL (-1 . 75,-1.75, . 125,DA,0. ,8) 

CALL  LINE(GAM2, FRQ, 200, 1,0,3) 


IN  Z-DIRECTION' ,-25,9.0,0.0, 


GAM 1(201} -GAM2 (201) 
GAM3(201)-GAM2(201) 

GAM1 (202)-GAM2(202 j 
GAM3 ( 202 ) -GAM2 ( 202 ) 

CALL  LINE(GAM1, FRQ, 200, 1,0,3) 
CALL  LINE (GAM3, FRQ, 200, 1,0,3) 


DO  55  K-l.NMODE 
DO  60  1-1,200 
CAM2(I)-TRAP(K,I) 

CONTINUE 

CALL  NEWPEN(2)  THIS  STATEMENT  CURRENTLY  FOR  COMMENT  ONLY 
CALL  LINE(GAM2, FRQ, 200, 1,0,3) 

CONTINUE 

ESTABLISH  ARRAYS  AND  RECORD  DATA  ON  PLOT. 


C(1 

C(2 

C(3 


-Cl 

-C2 

-C3 


Figure  16.  (continued) 
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202 


C 

C 

C 

C 

C 

8000 

8001 

8002 

9000 


9001 

9002 


9004 

9005 

9006 


9009 

9020 


C 

C 

C 

C 

C 

C 


R<l)-R01 

R(2)-R02 

R(3)-R03 

YD-6.5 

DO  202  1-1,3 

XD-8.0 

FLT-C(l) 

CALL  SYMBOL(XD,YD,0.125,'C-',0.0,2) 

CALL  NUMBER<999. 0,999. 0,0. 125, FLT, 0.0,1) 

FLT-R(I) 

XD-XD+1 .5 

CALL  SYMBOL (XD ,YD ,0. 125, 'RHO-' ,0.0,4) 

CALL  NUMBER(999. 0,999. 0,0. 125, FLT, 0.0,1) 

YD-YD-0 . 25 
CONTINUE 
XD-8.0 
YD-YD-O. 25 

CALL  SYMBOL(XD,YD,0. 125,  THICKNESS (M.)-  '0.0,15) 
CALL  NUMBER(999. 0,999. 0,0. 125, THICK, 0.0, 2) 

CALL  FINISH 
CALL  ENDOUT 


END  OF  GRAPHICS  COMMANDS 


FORMAT(F5.3,I2,I3) 
FORMAT(5E10. 3) 

FORMAT (3E1 0.3) 

FORMAT ( '  C1(M/SEC)-',E11.4,' 
'  *  THICK (M)-'  ,E10. 3) 


CE11.4, 


C2(M/SEC)=' , El 1.4,'  C3(M/SEC)> 


EPS-' ,F5. 3, '  NSIG-'  ,1 


FORMAT (  FREQUENCY  RANGE  IN  KHZ:  0-', El 0.3) 

FORMAT (  CONVERGENCE  CRITERIA  FOR  ZBRENT: 

C2 , '  MAXFN-'  13) 

FORMAT  C  MOdE  M-',I3,'  CUTOFF  FREQ-' , Ell .4) 

FORMAT ( '  R01-'  Ell. 4,'  R02-' ,E11 .4, '  R03-',E11.4) 

FORMAT C '  ZBRENT  FAILED  TO  CONVERGE  IN  MAXFN  ITERATIONS.  MAXFN-',I3 
C  '  I-'  15) 

f ORMAT ( '  ’  TOTAL  NUMBER  OF  MODES  M-',I3) 

FORMAT ('  ERROR . EXCEEDED  DIMENSION  FOR  TOTAL  NUMBER  OF  MODES.') 

STOP 

END 


EXTERNAL  FUNCTION  'F(X)'  IS  USED  BY  ZBRENT  TO  FIND  A  ROOT 
'B'  SUCH  THAT  'F(B)«0' . 


FUNCTION  F(X) 

REAL*8  X  F 

COMMON / AREA 1 /GAM II , GAM2I , GAM3I , RO l , R02 , R03 , I , THICK 

DU1-X**2-CAM1I**2 

DU2-GAM2I**2-X**2 

DU3-X**2-GAM3I**2 

IF(DU1  .LT.  0.)  WRITE(6,4000)  X,I,GAM1I 
IF(DU2  .LT.  0.)  WRITE(6 , 4001 )  X,I,GAM2I 
IF(DU3  .LT.  0.)  WRITE (6, 4002)  X,I,GAM3I 
ALPH 1-SQRT (DU  1 ) 


Figure  16.  (continued) 
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ALPH2«SQRT(DU2) 

ALPH3-SQRT(DU3) 

G-R02*R03*ALPH2*ALPH1+R02*R01*ALPH2*ALPH3 
H-RO 1*R03*ALPH2**2-R02**2 *ALPH 1 *ALPH3 

FORMAT T 
FORMAT (* 

RETURN 
END 

//DATA. INPUT  DD  * 

0.0  5200 

8.0  E03  0.5  E03  8.0  E03 

5.0E03  2.0E03  5.0E03 

/* 


4000 

4001 

4002 


DU2<0. 

DU3<0. 


X»'*E11.4/ 

x-*!eii.4/ 


GAM1I-*,E11.4) 

GAM2I-*,E11.4) 

GAM3I-*,E11.4) 


0.1E0  10.1E03 


* 


Figure  16.  (continued) 
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